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The mind of man plans his ways, 

But the Lord directs his step. 

Proverbs 16:9 
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Introduction 


Background 

The theory of aerodynamic sound was developed by Lighthill [ 1 ] through the introduc- 
tion of the acoustic analogy. This analogy views the limited fluctuating fluid region as an 
external stress field acting on the uniform acoustic medium at rest. Thus by rearrangement 
of the usual statements of mass and momentum conservation into an inhomogeneous wave 
equation, Lighthill found a general procedure for estimating the intensity of the sound 
produced in terms of the fluid flow. The acoustic analogy proceeds with the assumption 
that the localized flow field is unaffected by the acoustic field, completing the separation 
between aerodynamics and acoustics. While others have described the acoustic radiation 
as driven by distributions of sources, the acoustic analogy gave a consistent and rational 
way to determine the strength of these sources. Even though the analogy separates the 
acoustics from the full fluid dynamics problem, physical insight is needed in determining 
the necessary information from the aerodynamic flow. 

The extension of Lighthill’s theory to include surfaces in the flow was first made by 
Curie [ 2 ] in 1955, and generalized to sound produced by surfaces in arbitrary motion by 
Ffowcs Williams and Hawkings [ 3 ] in 1969. The Ffowcs Williams and Hawkings (FW-H) 
equation has found great success, especially as a tool for estimating the sound generated 
by propellers and helicopter rotors. This success may be partially attributed to Ffowcs 
Williams and Hawkings giving severed forms for the source terms. Using the theory of 
generalized functions in the derivation, Ffowcs Willi ams and Hawkings choose a source 
description of monopole, dipole and quadrupole terms, which have come to be known as 
thickness, loading, and quadrupole terms respectively. Although this particular description 
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is not unique, it yields meaningful results when each the sources are considered individually. 
For example in Farassat’s [4] early work, which used the FW-H equation as its foundation, 
only the thickness noise was calculated for helicopter rotors because of the difficulty in 
determining the aerodynamic field. Later the contribution of the loading source was added 
as experimental and computational results became available. 

The volume source, or quadrupole, has been a subject of lively discussion and specu- 
lation in rotor noise prediction often because of the difficulty in assessing its importance 
when the flow is complicated. It has been argued by Farassat /4j that the quadrupole may 
be neglected when the turbulent flow region is small, as in the case of very thin propellers 
or low speed rotors. However, probably the most fundamental reason that the volume term 
has not been included in many analyses is because it requires a detailed knowledge of the 
flow field around the body in advance. Until recently this information was nearly impossi- 
ble to attain. Indeed, many numerical prediction codes based on the FW-H equation have 
neglected the quadrupole source term precisely because of the inability to determine the 
correct source strength together with the analytical and numerical difficulties in dealing 
with a volume distribution of sources. 

As propellers and helicopter rotor speeds approach sonic speed, it is no longer prudent 
to neglect the quadrupole source, since any argument based on the normal ordering of 
monopole, dipole and quadrupole becomes meaningless. Hanson and Fink [5] were the 
first to demonstrate through numerical calculations the importance of the quadrupole 
source for high speed propellers, while Schmitz and Yu [6] have found similar results for 
helicopter rotors. Ffowcs W illiams [7] has shown that in some cases the quadrupole term 
may be the most significant source term. 
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Today prediction methods based on the FW-H equation are quite advanced. One good 
example is that of Schultz and Splettstoesser [ 8 ] where model scale helicopter rotor noise 
is compared to an approximate prediction of all three source terms of the FW-H equation. 
The aerodynamic blade pressure was measured on the blade and used as input in the acous- 
tic theory. Predictions were good for all the cases presented, which included high speed 
impulsive and blade vortex interaction noise. Unfortunately, when a prediction is needed 
we do not generally have the luxury of having measured blade surface pressures. Even so 
accurate results axe still necessary. Determining the time dependent surface pressure, not 
to mention the pressure and velocities of the fluid particles near the body, is difficult. These 
calculations are especially troublesome when the volume source is most needed; i.e. when 
the flow is transonic and sophisticated aerodynamic analyses are needed. Although current 
high speed computers and present day computational fluid dynamics (CFD) methods may 
allow one to predict the aerodynamic field, much of the advantage of formulating the acous- 
tic problem separately from the aerodynamics problem is lost. Indeed the development 
and operation of the CFD codes is a major task in its own right. 

Farassat [ 9 ] has proposed to reduce the quadrupole from a volume source to an ef- 
fective surface source by considering the time dependent shock surfaces as the primary 
contributor in the full volume source. Recognizing that the shock is a discontinuity in 
the fluid and guided by the mathematics, Farassat has shown that a shock surface noise 
contribution should improve the complete rotor noise calculation when shocks exist. While 
surface integration is preferable to volume integration in the acoustic calculations, this path 
demands more information from the aerodynamicist. In theory the entire flow field is not 
needed, but in practice the calculation of the full time dependent, three dimensioned flow 
field will only be the first step in determining the shock location, geometry, and strength. 
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Two alternate approaches to the FW-H for predicting the sound of moving bodies have 
been pursued in the past few years. The first alternative is to consider the FW-H equation 
as an integral equation and solve for the entire flow field. Long [10], Farassat and Myers 
[ll], Brandao/J 2,i3j , and Lee [14] have all attempted to solve singular integral equa- 
tions based on the FW-H equation with varying degrees of success. Hanson [15] and Das 
[16] have developed integral methods based on aeroacoustic ideas. Morino and his col- 
leagues/17,18/ have used a similar mathematical approach for propeller and helicopter 
aerodynamics. The acoustic integral equation has the advantage that the method used 
for aerodynamic and acoustic calculations is unified, yet it is apparently difficult to im- 
plement a nonlinear method along these lines. While this approach does not seem to take 
any advantage of recent progress made in CFD, it is hoped that much of the experience 
learned in developing noise prediction codes can be transferred directly to the develop- 
ment of aerodynamic codes based on the FW-H equation. It is also hoped that new codes 
developed along this line would not be significantly more complex than present day noise 
prediction codes. Ffowcs Williams [19] describes how this idea has been the motivation for 
acoustic diffraction problems, however, in general, diffraction theory has only been applied 
to simple geometries. 

The second alternative to the FW-H equation recognizes that a major effort is required 
to develop and use CFD codes just to determine the input for the FW-H equation. Since 
the acoustic fluctuations are part of the complete fluid dynamics problem, there is no 
fundament al reason not to solve the combined acoustics and aerodynamics field using an 
unsteady CFD code. This approach has the advantage that the nonlinear features in the 
flow, like shocks, are accounted for automatically along with local sound speed variation 
and convection effects. These effects are included in the FW-H quadrupole source, albeit 
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in a somewhat unnatural way. The difficulty with this approach is that non-physical 
boundary conditions and numerical dissipation are required for stable numerical methods, 
even though they may cause potentially insurmountable problems in the simultaneous 
calculation of aerodynamics and acoustics. Also the acoustic quantities are several orders 
of magnitude smaller than the aerodynamic quantities, a condition which implies that 
very accurate solutions axe required. To achieve a high level of accuracy, the flow field 
must be resolved sufficiently by the computational grid. Finally, only near field solutions 
will be possible in the foreseeable future since execution times and memory requirements 
inherently limit job sizes. 

Korkan, Von Lavante, and Bober [20] have overcome the near field limitation by com- 
bining a non-linear Euler calculation for a high speed propeller with a Kirchhoff formulation 
to determine the sound outside the computational grid. Thus linear acoustics is used in 
the region where appropriate and the near field calculation is completely nonlinear; all in 
the spirit of the acoustic analogy. Ffowcs Williams and Hawkings suggested that their 
formulation could be used as a Kirchhoff formula, but it took the advancement of both the 
digital computer and computational fluid dynamics for this to become feasible. 

The results of Korkan et. a l.[20] were not as good as one might hope, but their work 
did bring to the surface many of the issues, such as grid resolution, artificial viscosity, etc., 
which must be addressed. Purcell, Strawn, and Yu [21] as well as Lyrintzis and George 
[22] have applied the Kirchhoff approach to idealized helicopter rotor problems but there 
is some question about the extension of their work to complete three dimensional rotating 
blades. For this reason, Farassat and Myers [23] have given a new derivation of Morgans’ 
[24] original Kirchhoff formulation using the modern mathematics of generalized functions 
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and differential geometry, extending the result to include a control volume in arbitrary 
motion. . 

The numerical coupling of a nonlinear inner solution with a linear outer solution 
for combined near and far field acoustics prediction appears to be an efficient use of the 
numerical tools at hand. .This methodology depends crucially upon an inner solution which 
is accurate in both space and time within the entire control volume. The pressure, density, 
and fluid velocities along with their time and spatial derivatives at the control surface 
will be required as input into the Kirchhoff formulation. Cleaxly if there are significant 
damping or dispersion errors in a CFD solution, the outer acoustic solution cannot be 
accurate. Hence the coupled solution depends significantly on the inner solution. 

Outline of the Dissertation 

Two groups of problems are considered in this dissertation in order that some of the 
relevant issues in predicting the acoustic field of bodies moving at high speed may be 
examined. Since the quadrupole source has been shown to be important for accurate noise 
prediction but is both analytically and numerically difficult to evaluate, an analysis of this 
source term is made in Part A for the incompressible limit. While it is clear from Hanson 
and Fink. [5] that the quadrupole will be most effective when the flow field is transonic, it 
is difficult to calculate exact solutions for transonic flow. For the idealized case of inviscid, 
incompressible flow, however, the exact solution for a circular cylinder is well known and 
it is shown that the contribution to the pressure field by each of the thickness, loading and 
quadrupole terms can be calculated separately and exactly. In this case, it is shown that 
the lift on the cylinder due to circulation is divided equally between the loading source 
and the quadrupole source. 
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These results are generalized for the case of a Joukowski airfoil. Here it is found 
that that the volume source can be eliminated for incompressible flow by manipulating the 
quadrupole and changing variables. This new incompressible version of the FW-H equation 
is compared with thin airfoil theory showing the role of the quadrupole in an aerodynamic 
theory is to enable lift generation. The case of a vortex passing a circular cylinder is also 
considered, where it is found that the quadrupole is required to account for the concen- 
trated vorticity in the field. Keep in mind that in Part A the emphasis is upon the role of 
the FW-H source terms rather than the solution of the inviscid, incompressible flow past 
a cylinder. 

For the second part of this work, it is postulated that CFD will play an important 
role in the future of acoustics; either to provide input into the quadrupole source term 
or to calculate directly the near field acoustics. In Part B, a finite volume, multistage 
time stepping, Euler code is used to investigate the use of present day CFD algorithms 
for the direct calculation of acoustics. The two dimensional, compressible, inviscid flow 
about an accelerating or decelerating circular cylinder is used as a model problem. The 
time evolution of the energy transfer from the cylinder to the fluid, as the cylinder is 
moved from rest to some nontrivial velocity, is clearly seen. Energy is the main quantity of 
interest in the calculations since various components of energy have physical meaning. By 
examining the temporal and spatial characteristics of the numerical solution, a distinction 
can be made between the propagating acoustic energy, the convecting energy associated 
with the entropy change in the fluid, and the energy contained in the local aerodynamic 
field. Systematic variation of the cylinder acceleration shows that the radiated acoustic 
energy depends strongly upon the rate of acceleration or deceleration. 
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Longhorn [25] has found that for an impulsively started sphere there is an equipartition 
between the kinetic energy in the local flow field and the propagating acoustic energy, in 
the low Mach number limit. In our calculations for the impulsively started cylinder, we 
find that an equipartition between of the energy in local aerodynamic field and the energy 
transported to the far field through propagation and convection exists. This is an extension 
of the low Mach number results. When the Mach number exceeds the critical Mach number 
and shocks persist in the steady state flow field, there is no longer an equipartition of energy 
since the cylinder continues to do work on the fluid, offsetting the drag associated with the 
shock. The additional energy input into the fluid convects downstream from the cylinder 
along with the newly generated vorticity. 

Unlike the low Mach number situation, the energy of the aerodynamic field has a 
potential energy component, due to the significant second order compression of the fluid, 
in addition to kinetic energy. In our calculations, energy which leaves the control volume 
is divided between acoustic energy and energy associated with the change of entropy in 
the fluid. The computational grid has a large effect on the ratio of acoustic energy to 
entropy associated energy, while the role of the artificial viscosity seems to be of second 
order. We have found that the entropy terms were nearly negligible in all cases where 
the cylinder was started slowly. The calculations in Part B both give insight in to the 
interesting nonlinear interactions which occur when the cylinder is rapidly started as well 
as the numerical aspects of predicting the acoustics as part of the unsteady aerodynamic 
field. 
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Part A: Exact Calculation of Quadrupole Field 


Introduction 

In an effort to gain a new understanding about the FW-H quadrupole in both acous- 
tic and aerodynamic applications, we have chosen some simple problems for which the 
incompressible flow field can be determined analytically using the two dimensional veloc- 
ity potential. For a circular cylinder, each of the source terms are calculated separately and 
compared with the exact potential solution. The forces on the cylinder due to pressure are 
compared as well. This problem helps to explain the results and difficulties of Brandao’s 
application [12,13] of the FW-H equation to aerodynamics. 

The circular cylinder solution suggests a new description of the quadrupole term 
which is useful in identifying the volume and surface terms immediately from the exact 
solution. As an example, the problem of a circular cylinder moving near a vortex filament 
is examined to demonstrate the role of vorticity in the aerodynamic field. Following this, 
the contribution of the volume and surface sources are determined for a Joukowski airfoil. 
Each of these cases illustrate the various contributions of each of the source terms for 
incompressible flows. Finally the aerodynamic role of the quadrupole source is discussed 
along with a summary of Part A. Another analytical consideration of the role of quadrupole 
sources for exact compressible flow problems is given by Ffowcs Williams[7j. 
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The Incompressible FW— H equation 

The differential form of the Ffowcs Williams - Hawkings equation may be written as 


d 2 

dt 2 


, a 2 

’ dxidxi 


(p'HU)) = !(***/» - £<*«,*(/» + ^ 


(TijH(f)) (1) 


where c is the sound speed, p is the density, p' is the density perturbation, and the subscript 
o refers to the value of the variables in the undisturbed medium. Lighthill’s stress tensor, 
T{j, is defined 

Tij = puiUj + pij - c 2 p'Sij (2) 


where U{ is the fluid velocity, p{j is the compressive stress tensor, and is the Kroeneker 
delta. In equation (1) the derivatives are to be interpreted as generalized derivatives, while 
H(f ) and 6(f) are the Heaviside and Dirac delta functions respectively. The surface of 
the body is described by the equation / = 0 with / being defined such that V/ = n, the 
outward unit normal vector. Note that the monopole and dipole terms, commonly referred 
to as thickness and loading sources, act only on the surface / = 0 while the quadrupole 
source acts throughout the volume / > 0. The velocity term v n in the thickness source is 
the component of the body velocity normal to the body surface. 

For an incompressible, in viscid fluid, p' = 0 and p,j = p8{j . The FW-H equation may 
then be written 


VV#(/)) = -§ i (pvnS(f))+-~:(p'n i S(f)) - g^-(pu (3) 

The subscript o on the variable p has been dropped since the density is constant and the 
perturbation pressure p' is used rather than p to emphasize the fluctuation in pressure and 
simplify accounting later. This is allowable since V 2 p ; = V 2 p. Equation (3) can also be 
obtained by first using the isentropic relation p 1 = c 2 p' and then allowing the speed of 


10 



sound c to approach infinity to account for incompressibility. To obtain a solution to this 
equation, we need to determine the strength of each source term. 


Velocity Potential Solution for a Circular Cylinder 

The velocity potential for a circular cylinder of radius a, in a frame of reference in 
which the cylinder is moving, is unsteady and known to be 

„2 v-a 

= ( 4 ) 


Here r ,8 are the polar coordinates of x with respect to the cylinder center x 0 , v(<) is the 
cylinder velocity, K is the bound circulation on the cylinder, and r is a unit vector in the 
direction from the cylinder center to x. See figure 1. The pressure perturbation, given by 
the Bernoulli equation, is 


, 1 2 

P=P-P0 = --PU -P-Qi 


(5) 


where 


2 _ 


u — 


4 

* 2 


Ka 2 


vt + 


_K 2 _ 

47r 2 r 2 

\juf CL 2 f o 9-, (X 2 dv Kv* 


= -v*+ 3 

r* -kt* 


dt 


2 irr 


( 6 ) 

(7) 


Here v n = v • n and vt = v • t are the components of velocity normal and tangential to the 
surface, respectively. Note that for the circular cylinder r and n are equivalent while n and 
t are both unit vectors. Also notice that p 1 — ►Oasr— >ooas required by the physics of 
the problem. These terms are written out so that they may be compared with the solution 
gained from the FW-H equation. 
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f = 0 

Figure 1. This figure shows the geometry of the circular cylinder and the definitions of the 
variables used in the calculations. 

Acoustic Solution for a Circular Cylinder 

The solution of the incompressible FW-H equation (3), which we will call the acoustic 
solution, can be obtained using the two dimensional free space Green’s function of the 
Laplace equation, G(x; y) = In |x — y |/27 t, since the pressure and velocity fields are known 
and the geometry is simple. The Green’s function solution for the model equation 

V 2 <£(x,t) = Q(x,i) (8) 

is given by 

oo oo 

2n<f>(x.,t) = J j Q(y,t)\n\x-y\dy. (9) 

— OO — OO 

The free space Green’s function is appropriate since the source Q(x,t) is a generalized 
function which extends the ordinary function to the entire two dimensional space. The 
presence of the generalized functions 6(f) and H(f) remind us that the source terms in 
the FW-H equation have been extended as part of its derivation. Faxassat discusses the 
powerful technique of extending ordinary functions in references [4,26] . 
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In Appendix Al, the integration for each of the source terms in the FW-H equation 
is carried out analytically. When this is done, the pressure perturbations obtained are 


and 



( 10 ) 

( 11 ) 


. P [<T ? Ka? 

^ = - 21 ^ + ^ V< + 


J£!_i 

47 r 2 r 2 J 


( 12 ) 


Here the subscripts t,l, and q on p' refer to the thickness, loading and quadrupole con- 
tributions, respectively. The subscripts n and t on v refer to the normal and tangential 
components of the cylinder velocity. It is immediately clear when comparing these equa- 
tions with the potential solution, equations (5-7), that the thickness and loading sources 
correspond exactly to —pd<f>/dt and the quadrupole contribution corresponds to —pu 2 /2. 
This finding warrants further exploration to determine if this correspondence will be true 
in general or if it depends upon the unique geometry of the circular cylinder. Remember 
that these results apply when (f> is defined for the case of a moving body in a stationary 
fluid. 

Notice that the complete fax field solution is given by the thickness and loading terms, 
but if there is circulation, K ^ 0, and no acceleration of the cylinder, dv/dt — 0, the 
quadrupole contribution can be as important as the thickness term. The quadrupole serves 
to provide a near field pressure adjustment to the thickness and loading pressures. In effect 
the linearized pressure perturbation, p' = —pd<j>/dt, is determined by the thickness and 
loading sources alone, while the quadrupole accounts for the local kinetic energy in the 
fluid. Figure 2 shows the relative contributions of each of the source terms for a cylinder 
with circulation. This incompressible solution, which contains no acoustics, highlights that 
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it is possible to use the FW-H equation to predict the aerodynamic pressure field as well 
as acoustic pressure perturbations. 


Force on the Cylinder 

The force on the cylinder can now be calculated by integrating the pressure over the 
cylinder surface. The force per unit length is found to be 


F = F* + Fj + F 9 


where 


Ft 
F / 
F 


1 dv 
— -m— 

2 dt 

1 dv 1 f 

~2 m Tt + r K{y * k) 


? — 


\p R iy x k) 


( 13 ) 


( 14 ) 

( 15 ) 

( 16 ) 


Here m — pna ? which is the added mass of the cylinder and k == n x t. The force 
composed of Ft and the first term of Fj is an inertial force needed to accelerate the added 
mass while the force composed of F g and the second part of F i is a lift force due to 
circulation. The inertial force is independent of the quadrupole, due to the symmetry of 
the cylinder, but one half of the lift force is given by the quadrupole term. This implies 
that if the FW-H equation is to be used for aerodynamic calculations, the quadrupole may 
be important for lifting problems. 
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c) Loading contribution d) Quadruple contribution 

Figure 1. The pressure perturbation around a circular cylinder, radius a — 1.0, moving at 
a velocity v = 1.0, and with circulation K = 7r. 
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A New Quadrupole Description 

Before a more definitive statement is made, let us first return to examine the way in 
which the quadrupole term was simply related to pu 2 / 2. In Appendix A2 it is shown that 
with no loss of generality the quadrupole term in the incompressible FW-H equation can 
be rewritten 

fod =V 2 {ipu J /f(/)} +/>v ■ {< X u H(f)} 

+ pV{(» n u-i« 2 n)«(/)} (17) 

where ( = V x u is the local vorticity of the fluid and p is constant. The surface term 
arises from the gradient of the generalized function H(f). The second term on the right 
hand side is zero for an irrotational flow ( V x u = 0). This form of the quadrupole clearly 
shows that if we only consider the linear solution, i.e. drop the pu 2 / 2 terms, that we still 
have a volume source and a surface source remaining. Hence the linearized solution is not 
equivalent to neglecting the quadrupole. This is not clear upon examining the original 
form of the quadrupole. The purpose of deriving this quadrupole expression is to separate 
explicitly the pu 2 / 2 part from the other parts, but (17) does have the advantage of showing 
the role of vorticity in the aerodynamic field. 

It is also useful to rewrite the thickness term 

£ W(/)} = -V • {pv n v6(f)} + p ^ n 6(f) (18) 

which is a form similar to the previous equation. The derivation of this relation is found in 
Appendix A2 as well. The incompressible form of the FW-H equation may now be written 

V 2 {p'tf(/)} = - V 2 {i pu 2 H(f)} - ,>V{(< x u)ff(/)} - V . {pv n ( u - v)S(f)} 

+ V •{(?' + i^ 2 )n S(f)}- P ^ n6(f). (19) 
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This new equation displays elements similar to Powell’s theory of vortex sound [27] where 
the quadrupole source region is identified with the vorticity of compact eddies in the flow. 
The second source term in equation (19) is restricted to the region in the flow where the 
vorticity is nonzero, while the third source term is written in such a way that it may be 
interpreted as a vortex sheet of strength t it — vt over the surface, since u — v = (ut — v*)t 
on the surface. For the important case of irrotational flow the volume source of this new 
equation is exactly —pu 2 / 2, as was the case for the circular cylinder. In the following 
problems it will be possible to calculate the exact potential solution and then directly 
identify the volume and surface components in the modified FW-H equation, equation 
(19). 

It is well to point out that if the body rotates this will affect both the first and last 
surface terms since the velocity v will have an additional rotational component. This 
explains the apparent discrepancy for the case of a circular cylinder rotating about its axis 
while also translating. If the cylinder rotates, the first surface source term appears to have 
a strength which depends upon the rotation rate of the cylinder, but the change in the last 
surface term exactly accounts for this apparent ambiguity. This must be the case for an 
inviscid fluid. 

Upon examination of the modified FW-H equation, we see that a change of vari- 
ables from p 1 to B = p' 4- pu 2 / 2 = —pd<f>/dt is appropriate. This choice of variables, 
which is related to the variable Howe [28] used for his nonlinear analogy, transforms the 
FW-H equation into 

V 2 (£ff(/)) =V • { (BA - pv„(u, - »,){)<(/)} - P^ ■ *«(/) 

-pV.((Cxu)ff(/)). (20) 
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Here the variable B is actually the linearized pressure perturbation. The linear contribution 
of the original quadrupole is now represented by the vorticity term and the pv n ut surface 
term. Notice that for an irrotational flow, the volume source term disappears, a reasonable 
result since the Laplace equation can be solved uniquely from the boundary data. 

An alternate derivation given in Appendix A2 makes it easy to see that we can also 
write 

V 2 (Bff(/)) = V . (BnS(f)) - p ^ ■ n S(f) - pV ((< x u )#(/)). (21) 

If we now recognize that V<j> = u and B = —pd<f>/dt, we find that (21) can be written in 
the form 

V 2 (BB(/)) = V . (BAS(f)) + VB . ftf (/) - pV ■ ((< x u)B(/)) . (22) 

This version shows that when there is no vorticity in the fluid, the solution can be expressed 
in terms of B on the surface alone, but it does not give much physical insight into the 
solution. This equation can be manipulated to give similar forms which are in terms of 
the fluid vorticity and B on the surface. 

Solution for Circular Cylinder with a Vortex 

Now we will consider a circular cylinder moving past a vortex to demonstrate the 
utility of arranging the source terms as a combination of surface terms, kinetic energy 
term, and vorticity term in equation (19). This example is chosen to demonstrate the 
role of the volume source terms when a region of vorticity exists in the aerodynamic field. 
For this problem, the complex velocity potential can be found using the Milne-Thomson 
circle theorem [29] to be the sum of the velocity potential of the cylinder alone, the free 
vortex alone, and an image vortex system. If z\ is the position of the free vortex expressed 
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as a complex variable, z\ = x\ + iyi, the image vortex system consists of a vortex of 
equal strength to the free vortex at the cylinder center, and a vortex of equal strength 
and opposite sense at the image point z 2 = a 2 /zi. Here the overbar indicates the complex 
conjugate of the complex number. The complex velocity potential for the problem is 
written 

= zXfi. + »(^ + r) ln (*) , trin(z-z!) _ tITn(z - z 2 ) 3 

z 27r 27 t 2tt 

after dropping a constant. In this expression the cylinder has a velocity of magnitude v 

directed at an angle a to the x axis, expressed as the complex velocity V = u(cos a+t sin a). 

The pressure perturbation p' is, from Bernoulli’s equation, equal to —pd<}>/dt — pu 2 /2, i.e. 

, 1 ,dw, 2 / dw . 

” = “2' , l9rl -' ,Re( ar ) (24) 

where 


6w Vo? i(K + T) 

rt O 1 n I 


2i rz 27r(z — zi) 27 r(z — z 2 ) 


^ _ _v— - —( Yl Y 2 ^ (oa\ 

dt~ dz 2t t K (z-zi) ( z-z 2 ) } ‘ K J 

Here dw/dz is the complex conjugate of the fluid velocity written in complex form, i.e. 
dw/dz = u\ — in 2 , and Vj, V 2 are the complex velocities of the free and image vortices. 
Although it is convenient for our purpose to write the velocities of the vortices as separate 
quantities, they are not independent. The velocity of the free vortex is that which a fluid 
particle at the position of the vortex center would have due to the combined velocity field 
of the cylinder and image vortex system, while the velocity of the image vortices is related 
to both the position and velocity of the free vortex relative to the cylinder. 
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This particular problem highlights a situation where the vorticity term in the volume 
source must not be neglected. Since the vorticity is concentrated at the point xi, the 
vorticity vector £ can be written -r<5(x - Xi)k where T is the strength of the vortex and 
k = n x t is a unit vector in the z direction. The pressure contribution of the free vortex 
may then be written 

p' v = ^-V- l pTk x uS(y -xi)ln|x-y|dy = prk x u| • - Xl ^ (27) 

Z7r J xi 2 tt x — xi r 

/> 0 

which may be rewritten in terms of complex variables as 

?; - = ' Re <2^rr^)>' (28) 

This can be recognized immediately as part of d<f>/di in equations (24) and (26). For 
this problem, the linear pressure perturbation, -pd(f>/dt, is the sum of the surface terms 
and the free vortex alone. The surface source contribution differs from the cylinder alone 
solution by exactly the contribution of the image system of the free vortex. Powell [30] 
and Ffowcs Williams [31] have shown this is true for turbulent boundary layers on a plane 
boundary as well. Notice that the contribution from the vortex alone in the flow field 
is part of d<f>/dt rather than pu 2 / 2. Thus the role of the volume source is primarily the 
inclusion of the free vortex pressure perturbation and secondarily to provide a nonlinear 
coupling of the various source components through the kinetic energy of the fluid. This 
suggests approximating the solution by the linear pressure perturbation, which neglects 
the nonlinear pu 2 / 2 part of p' . Figure 3 shows this approximation compared to the exact 
solution. The figure shows that this approximation varies little from the complete solution 
except in the vicinity of the free vortex and near the cylinder surface. 
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Figure 3. Comparison of p 1 and —pd^/dt for a circular cylinder moving in proximity to a 
free vortex. ( V = 1.0, K = 2.0, T = —2.0, Vortex position (ri,0i) = (2., 15 deg.) ) 



Solution for the Joukowski Airfoil 

Now as a further example consider the case of a Joukowski airfoil in an incompressible 
flow. The exact solution is readily obtained using the Joukowski transformation, ( = 
z + 1/z, to transform the complex velocity potential w(z) for the circular cylinder. The 
pressure perturbation p 1 can be written 

P' = -\p\^\ % - PMvty W 

where the airfoil has a steady complex velocity V and dw/d( is the conjugate of the 
complex fluid velocity. The value for dw/d( is obtained by transforming the solution for 
the circular cylinder with a free stream moving past into the £ plane and then subtracting 
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the freestream velocity. The freestream is removed since the velocity potential must be 
given in a frame of reference where the body is moving in a still fluid. Once this is done 
we find 


dw V a 2 iK 

d( (z - z 0 ) 2 + 2t t(z - z 0 ) 


- ) /— ~v 

- z„\' I dz 


(30) 


where V is the conjugate of V and z 0 is the center of the circular cylinder. Clearly an 
analytical term by term integration of the original FW-H equation using the potential 
solution is difficult for even this simple airfoil geometry. There is no vorticity in the flow 
field, therefore we can use the new quadrupole description to identify the volume and 
surface contribution to p 1 immediately from Bernoulli’s equation. Unfortunately such a 
separation is not possible with the original form of the FW-H equation. 

In Figure 4 the surface source contribution is compared with the velocity potential 
solution for a cambered airfoil at angle of attack. In this figure, it is evident that the 
volume source, -pu 2 / 2, has a negligible effect away from the airfoil. In figure 5, the 
pressure coefficient on the surface of the airfoil is compared with the surface and volume 
components, —pd(f>/dt and —pu 2 /2. The volume source contribution is most important 
near the forward stagnation point, and we can see by comparing the difference in pressure 
on the upper and lower surface of the airfoil for each component that the volume source 
term plays a significant role in the lift generation. 


Aerodynamic Implications of the Quadrupole 

From what we have already seen from the three previous examples, it is clear that the 
volume source term may be safely considered to have negligible influence on the pressure 
perturbation far away from the body. Recall however that the FW-H quadrupole accounts 
for half of the lift on a circular cylinder with circulation and the volume source term in the 
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Figure 4. The pressure perturbation components for a flow (v = 1.0) about a Joukowski 
airfoil (o = 1.13, z 0 = —.11 + .10t) at a = 5deg. 


Joukowski airfoil problem produced a smaller, yet significant part of the total lift. These 
apparently contradictory observations indicate that the quadrupole is important in any 
aerodynamic theory based on the the aeroacoustic approach. Without going far from our 
central purpose it is possible to shed some fight on the role of the quadrupole source in 
aerodynamics. 

As we begin to consider aerodynamics, it is important to realize that the primary 
purpose of aerodynamics is to determine the loading on the body, rather than to determine 
the perturbations far away from the body. For incompressible flow especially, it is obvious 
that it would be more straightforward to solve the Laplace equation for either the velocity 
potential or the acceleration potential using Green’s theorem to obtain the aerodynamic 
solution. The attraction of an aeroacoustic based aerodynamics theory is that it would be 
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Figure 5. A comparison of pressure coefficient contribution from the surface source, Cp 3 , 
and the volume source, Cp v with the toted pressure coefficient, Cp = 2 p' /(pU^), for a 
Joukowski airfoil (a = 1.13, z 0 = —.11 + .10*) at a = 5deg. 


unified with the acoustics, and hence one theory and one computer code could be used to 
solve the entire aero/acoustic problem. It is also possible that in developing such a theory 
new understanding may be gained by using an acoustic point of view. Further, the theory 
may be more general in nature, more robust, or its implementation numerically superior 
in some way. 


24 


Long [10] and Farassat and Myers [11] were the first to use the FW-H equation for 
aerodynamic problems. Long started with a linearized derivation of the FW-H equation in 
which he eliminated the quadrupole from the start, as is common practice in the prediction 
of the acoustic field of rotating blades. He was able to develop a singular integral equation 
which was useful for the thickness problem, but Long was forced to resort to an ad hoc 
numerical conditioning to achieve a lifting solution. Curious about this result, Farassat 
and Myers/1 1 J have shown that the angle of attack problem becomes an eigenvalue problem 
in Long’s case and cannot be solved, even though the unsteady aerodynamics problem is 
well posed. This finding seems to be supported by our earlier incompressible result where 
we have found that quadrupole did not contribute to the force on the cylinder due to 
unsteady motion. Brandao[12,13j has also developed an aerodynamics theory based on the 
FW-H equation but has reduced most of his work to the case of an incompressible fluid for 
simplicity. Since there has been some confusion about the results of aerodynamic theory 
based on the FW-H equation, we provide some insight into these problems in this section. 

The steady incompressible form of Brandao’s [12] aerodynamic integral equation may 
be written 



p'h • f 


dS = 



pv n \ ■ 
7*2 



(31) 


where f = (x — y)/|x — y|, r = |x — y|, and dS is an element of the three dimensional 
surface / = 0. This is a singular Fredholm integral equation of the second kind for the 
unknown variable p'. Equation (31) is also the incompressible limit of Long’s equation. 
We can develop a similar integral equation based on equation (20) by using the three 
dimensioned Green’s function for the Laplace equation for unbounded space. Having done 
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this, the steady form becomes 


_ . » If Bn ■ r 1 f 

f = 0 /=0 


pv n (u - v) • r 


dS 


-- / 
47T J 

/> 0 


p(C X u) r 


rfV. (32) 


This equation is valid for steady motion in an incompressible fluid and is closely related to 
the acceleration potential method. It differs from the normal acceleration potential theory 
in that it also accounts for the wake of a lifting body through the volume integral on the 
right hand side. In Appendix A3 it is demonstrated that classical thin airfoil theory can 
be obtained from the two dimensional form of equation (32) for both thickness and lifting 
cases. With this knowledge, we can compare (32) with Brandao’s equation to determine 
where the quadrupole is important. 

Examining (31) and (32), we can see that if we use the linearized pressure, i.e. if we 
let B — > p', the two equations axe similar except for the absence in Brandao’s equation 
of the u • f term in the first integral on the right hand side and the lack of any volume 
integral. As we noticed earlier, this combination of u — v can be thought of as a vortex 
sheet around the body which, in the three dimensional case, trails off the body forming 
a wake. The strength of this sheet is responsible for the lift generation, a feature absent 
when the quadrupole is neglected. 

Brandao[13/ recognized that the induced velocity acting through the quadrupole must 
be important for the lifting problem and used this as justification for an alternate inter- 
pretation of v on the surface. Brandao replaced v • f with (v n n + qt) • r, where q was 
determined iteratively during the numerical solution process. Having done this, Brandao’s 
results improved somewhat for blunt bodies with no lift, but the results were still unsatis- 
factory for lifting bodies. Equation (32) should have no such difficulty since it is shown, in 
Appendix A3, to give Munk’s Vortex-sheet theory [32] result for lifting thin airfoils. For 
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the two dimensional case, the volume integral is not needed if the wake is idealized by the 
Kutta condition. 

While equation (32) has been useful in the present discussion, in practical application 
it has the disadvantage that the fluid velocity u on the surface is regarded as being given. 
In general u will not be known and cannot be measured on the surface since real fluids 
obey the no-slip condition on the surface of a body. This can be overcome if equation 
(22) is used as the starting point since it is written in terms of the variable B alone on 
the surface. Since this equation is equivalent to the acoustic equation used in Appendix 
A3, it should also reduce to slender body theory and Munk’s vortex sheet theory under 
the assumptions of thin airfoil theory, even though it is more general than either of these 
theories. A further drawback, however, is that this form is only valid for an incompressible 
flow and it is not clear whether a compressible analog can be derived with the desired 
utility. In addition, an acoustic approach is probably complicated for the transonic case 
when the acoustic and aerodynamic fields are no longer distinct. In transonic flows, finite 
strength shocks may generate vorticity and the local speed of sound will vary. Accounting 
for these important effects will require the use of the full volume source term. Thus it is 
probably more natural for transonic flows to consider the numerical solution to the Euler 
or Navier-Stokes equations directly. 

Conclusions 

The aim of part A has been to gain more understanding of the importance of the 
quadrupole source in the FW-H equation. By using the incompressible flow solution, 
the exact contribution from each of the source terms in the FW-H equation was seen for 
both thick and thin bodies. For the special geometry of the circular cylinder, the sum of 
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the thickness and loading contributions to p 1 is proportional to d<j>/dt, when the velocity 
potential <f> is written in a frame of reference fixed to the undisturbed medium. The 
quadrupole contribution is —pu 2 / 2 in this frame. 

In general, the quadrupole source of the FW-H equation cam be divided into a group 
of surface terms and volume terms. The particular form of the surface source terms given 
in the new description is similar to that of the thickness and loading source terms in the 
original FW-H equation. The new volume source terms are related to the vorticity and 
kinetic energy of the fluid. This rearrangement of the FW-H quadrupole clearly shows 
that the quadrupole source has a contribution to even the linear pressure perturbation 
through the surface source terms and the vorticity term. In the special case of the circular 
cylinder, only the nonlinear, kinetic energy term of the quadrupole is nonzero, however. 

From the incompressible problems considered we can conclude that the quadrupole 
source term in the FW-H equation has two primary functions. As the only volume source 
term, the quadrupole is obviously necessary to account for any vorticity in the fluid. This 
was evident for the vortex-cylinder interaction considered. The quadrupole also provides 
a nonlinear coupling through the kinetic energy term. This is a localized effect, which 
nonetheless is important when the FW-H equation is applied to the aerodynamics problem 
since it accounts for a signficant portion of the steady lift. 

It is significant that both slender body theory and vortex sheet theory can be derived 
from the new incompressible form of the FW-H equation because this shows that the 
full incompressible form of the FW-H equation may be considered as the fundamental 
governing equation. Further development would be required to solve the singular integral 
equation in Appendix 3, even though Long, Farassat and Myers, and Brandao have already 
laid the groundwork. Developing an analogous compressible source description would be 
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much more complicated. A consideration of the compressible flow field is necessary to gain 
further understanding, so in Part B a numerical study of the compressible flow about a 
circular cylinder is carried out. 
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Part B: Numerical Calculation of Sound 
Generated in a Compressible Fluid. 


Introduction 

While the exact analysis of an incompressible flow field is useful, many of the important 
fluid dynamics phenomena can only be studied in a fully compressible flow. Transonic 
flows in particular are difficult to calculate analytically with any generality, yet with the 
computational techniques and computer power available today we can learn a great deal 
about complex, nonlinear flow fields by using computational fluid dynamics (CFD). In 
part B we have numerically solved the Euler equations, together with the continuity and 
energy equations, for the two dimensional problem of a circular cylinder accelerating from 
rest and the related problem of a decelerating the cylinder. The Euler equations were 
chosen because we specifically want to consider the case of an ideal, inviscid fluid so 
that the artificially viscous nature of the discrete numerical algorithm may be observed 
independently. The Euler equations are also useful because they allow us to consider 
transonic flow fields containing shock waves. 

A brief discussion of the numerical procedure is followed by an investigation of an 
accelerating cylinder. The transient process of starting the cylinder impulsively to Mach 
0.4 is presented in detail, then other starting scenarios are considered. In this analysis, 
we find that the energy of the fluid is a useful quantity to track since we are ultimately 
interested in the acoustic energy flux generated by accelerating the body. As we consider 
energy, Myers’ exact energy corollary [ 33 ] suggests physically meaningful partitioning of 
energy. When the energy is separated into kinetic, potential, and entropy terms along 
with the flux of these terms out of a control volume, we are able to observe the transfer 
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of energy from one form to another as a function of time. When a cylinder is started 
impulsively, we find an equipartition between the local kinetic and potential energy and 
the flux of acoustic and entropy energy out of the near field. This equipartition of energy is 
evident until the velocity of the cylinder reaches the critical Mach number, at which time 
the flow field contains shocks in the steady flow field. When this happens, the cylinder 
must continually input energy into the system to offset the shock induced drag. 

The case of a cylinder decelerating to rest is also considered. We find that if the 
cylinder is stopped impulsively, nearly half of the kinetic energy is rapidly transformed into 
potential energy as the energy in the local field propagates away as sound. Unfortunately a 
significant portion of the acoustic energy propagating away from the cylinder is converted 
to entropy energy, if the rate of deceleration is high. This is probably a numerical limitation 
which can be identified for the high acceleration cases as well. 

Normally entropy energy, i.e. the energy associated with increases in the fluid entropy, 
is due either to an increase in entropy across shocks or to fluid viscosity. Since we have 
an inviscid fluid, only entropy generated by shock waves has physical meaning. In the 
fined section of part B, we consider the potential sources of the entropy energy in our 
calculations. We find that when the final velocity of the cylinder is near the critical Mach 
number, a shock forms during a rapid start up process and entropy is generated. We also 
demonstrate that the implicit damping - due to the discrete nature of the algorithm, and 
the explicit damping - added for numerical stability, are both numerical mechanisms for 
the transfer of energy into the entropy energy term. We have found that the implicit 
damping related to the grid resolution is of primary importance and that the entropy 
energy term gives a quantitative indication of the validity of the numerical solution for the 
sound. 
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Numerical Method 


The Euler equations were chosen as the appropriate set of governing equations since 
they represent the physical characteristics of both sound and shock waves. The Euler 
equations also allow rotational flows to be considered if vorticity is generated by a shock 
or input into the system as an initial condition. In addition, the artificial viscosity of 
the discrete numerical solution is suspected of adversely affecting the calculated acoustic 
solution, therefore it is best to treat the inviscid problem first. With the assumption of 
an inviscid fluid, we need to remember that the solution for the circular cylinder will not 
be physically realizable since viscous separation is fundamental in the real problem. The 
advantage, however, is that any viscous phenomenon observed in the calculations must be 
related to the discrete nature of the numerical solution. 

A finite volume method with multistage time stepping of the Jameson typ e[34,35] 
was programmed to solve the two dimensioned Euler equations. This method is a cell 
centered scheme which is second order accurate in space and time in this implementation. 
The central difference methods of this type produce good solutions in smooth regions of the 
flow, but axe prone to oscillations in the neighborhood of shock waves. These oscillations 
require the explicit addition of dissipative terms. In the present code, Jameson’s adaptive 
blending of second and fourth order dissipation is used. The fourth order term provides 
a low level of background dissipation while the second order term provides a higher level 
of dissipation near shocks. The adaptive nature of the dissipation comes through the use 
of a switch, based on the pressure gradient, which turns on the second order dissipation 
when the gradient becomes large. The action of this adaptive dissipation on propagating 
acoustic waves is not well understood. 
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A five stage, Runge-Kutta like, time stepping scheme is used with the dissipative 
terms evaluated only in the first two stages to reduce the number of computations. This 
combination is second order accurate in time. Non-reflecting boundary conditions, based 
on the extrapolation of the linearized characteristic variables, are used at the outer com- 
putational boundaries. Symmetry is assumed between the upper and lower half planes. 
Therefore computations are carried out only in the upper half plane and symmetry bound- 
ary conditions are applied on the x axis. The surface boundary conditions consist of a 
statement of no flow through the body surface together with a linear extrapolation of the 
pressure onto the surface. 

One slightly unusual feature of the Euler code developed for this work is that the 
computational grid is moving rather than fixed. Thus, unlike most two dimensional aero- 
dynamic codes, the body is moving through a stationary fluid and the dependent variables 
are specified in a reference frame fixed to the undisturbed medium. The Euler equations 
appropriately modified for a moving control volume are written in vectorial form as 




^ + V x -F = 0 


(1) 

where 

( ' 1 
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which is the conservative form of the Euler equations. Here U is the vector of dependent 
variables - density, momentum, and total energy; q is the fluid velocity (u, v); v is the 
velocity of the grid - normally the body velocity; and i, j are the unit vectors in the 
coordinate directions. In these equations, the subscript X indicates that the spatial inde- 
pendent variables are moving with the control volume. The UtT term in the flux matrix F 
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arises from the time derivative of the moving coordinates. Since these additional terms will 
always be apparent in the flux matrix when the control volume is moving, the subscript 
X will not be included in subsequent equations. A moving grid was chosen to simplify 
the treatment for an accelerating body and to cast the problem in a frame of reference 
familiar to acousticians. The development of the equations for a moving control volume is 
addressed in more detail in Appendix B. 

The calculations presented have been nondimensionalized; length by the cylinder ra- 
dius, density and sound speed by the their values in the undisturbed fluid. With this 
particular nondimensionalization the pressure in the undisturbed medium is I/ 7 , where 7 
is the ratio of specific heats, and the speed of the body is equivalent to the Mach number 
of the body. An advantage of this system is that sound waves travel at a velocity close 
to unity plus the local convection speed. Thus it is fairly easy to determine where wave 
fronts should be expected, especially for an impulsive motion of the body. 

The finite volume, multistage time stepping method of Jameson was chosen for this 
work based on several considerations. One attraction is that this method has been used 
in the CFD community for several years and is efficient and robust. Secondly, the various 
parts of the method are distinct from each other and therefore may be changed indepen- 
dently. For example, there are many alternate time stepping procedures available if one 
changes the number of stages and the weighting coefficients of each stage. Indeed the sta- 
bility and efficiency of the algorithm can be changed considerably while still keeping second 
order, or higher, time accuracy. Likewise, Jameson [36,37] has shown that the algorithm 

can be transformed into an upwind or total variation diminishing (TVD) formulation by 
changing the form of the explicitly added dissipation. In addition, the flux balance can 
be changed from the cell centered, central difference scheme to a nodal scheme and there 
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is no particular difficulty in using this method on unstructured grids. This flexibility has 
not been utilized in this work, but it was a driving factor in choosing this method and 
will benefit future research using this code. A more detailed description of the precise 
numerical implementation is given in Appendix B. 

The Accelerating Cylinder 

The process of accelerating a body can be a potent noise source if the rate of accel- 
eration is sufficiently large. It is natural then to consider the transient behavior of the 
compressible flow field in this situation. Several model problems have been studied analyt- 
ically for various limiting conditions. Ffowcs Williams and Lovely [ 38 ] considered the case 
of a sphere that is suddenly brought into a low Mach number translation in an inviscid, 
compressible fluid. G. I. Taylor [ 39 ] also considered a sphere which, after Jin impulsive 
start, decelerated due to the drag induced by the transient pressure field. In each case, 
there is an equipartition of energy between the kinetic energy in the local flow around the 
body and the radiated acoustic energy. Longhorn /25j also found that the work required to 
start a sphere impulsively was twice the amount needed if the sphere was started slowly. 
The additional work required to rapidly accelerate the sphere supplies the acoustic energy 
which radiates as an intense sound. Each of these analyses are based on low Mach number 
motion, a limiting assumption we wish to remove. 

Hill [ 40 ] carried out a similar analysis determining the linear wave field induced by the 
sudden starting of an infinite wavy plane. He also found an equitation of energy between 
the steady evanescent field following the surface and the outwardly traveling sound for the 
case of subsonic surface phase speed. Hill’s problem displayed rather interesting transient 
behavior for the cases of subsonic, supersonic, and exactly sonic phase speeds. Nevertheless 
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none of the analytical approaches are able to adequately include the important nonlinear 
effects found in transonic flow. 

Description of the Starting Process. 

With these examples in mind, we shall first consider the transient flow field around 
a circular cylinder which has impulsively started in Mach 0.4 leftward translation. By 
examining this case in some detail we will be able to interpret the results for other Mach 
numbers and accelerations. The computational grid used for these calculations is a polar 
grid which models half of the flow field. The azimuthal direction is divided into 95 cells 
while the radial dimension of the grid cells increases linearly up to ten cylinder radii from 
the center of the cylinder, after which the radial increment remains constant. Figure 
1 shows the inner part of the computational domain. Notice in figure 1 that the grid 
resolution is very fine, especially near the cylinder surface. 



-10 -5 0 5 10 


Figure 1. This is the inner part of the computational grid used for the numerical calcula- 
tions. The grid is divided into constant azimuthal increments of A 9 = .0331 radians and 
linearly increasing A r from the surface to 10 radii from the cylinder center, after which 
A r remains constant. The initial value of Ar = .0397 near the cylinder surface. The x 
axis is shown in the figure to provide a spatial scale of reference. 
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When the cylinder impulsively starts moving leftward, an expansion wave propagates 
behind the cylinder while a compression moves ahead of it. This is seen in the time history 
sequence of density perturbation contours plotted in figure 2. In figure 2 the density 
perturbation p 1 is scaled by r 1 / 2 , r being the distance from the cylinder center, to offset the 
effect of cylindrical spreading. The expansion moves away from the cylinder more quickly, 
due to the motion of the cylinder. The steady flow field around the cylinder is nearly 
established by t=20.0. A careful look at figure 2 will also reveal that the rightward moving 
expansion is followed by a weaker recompression and the leftward moving compression by a 
weak expansion. The calculations show evidence that the signal has a sharp wavefront and 
the characteristic tail of cylindrical waves. Note that in all the contour plots the maximum 
value on the legend does not indicate the absolute maximum in the plot, but the regions 
colored with the maximum color should be understood to be greater than or equal to the 
maximum on the legend. Likewise the minimum color should be interpreted as less than 
or equal to the legend minimum, although some effort has been made to ensure that the 
maximum and minimum values in the contour region are not greatly different from the 
legend range. The values plotted in the contour plots are at the center of the grid cells, 
hence in some of the figures a small wedge shaped void is noticeable along the x axis. 
This is an artifact of the plotting. Contour plots of the pressure perturbation p 1 are nearly 
indistinguishable from the density perturbation and hence are not plotted. 

The surface pressure coefficient is plotted in figure 3. In the figure, — C p is plotted 
against x so that the negative pressure peak is plotted upward, as is commonly done in 
aerodynamics. Immediately after the impulsive start the pressure on the surface corre- 
sponds to p' = p 0 c 0 Vn , which is expected for an impulsive motion. The transition from 
the asymmetric to the symmetric pressure distribution as the flow reaches steady state 
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Figure 2. This figures shows a time sequence of density perturbation p' contours for the 
cylinder impulsively started to Mach 0.4. The density perturbation has been scaled by 
r 1 / 2 to account for cylindrical spreading. Here r is the distance from the cylinder center. 


38 


ORIGINAL FAGE IS 
OF POO** QUALITY 











k) t-20.0 


I) t-22.0 


Figure 2. (continued) This figures shows a time sequence of density perturbation p ' contours 
for the cylinder impulsively started to Mach 0.4. The density perturbation has been scaled 
by r 1 / 2 to account for cylindrical spreading. Here r is the distance from the cylinder center. 
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is seen in figures 3a-l. In the density perturbation contours it was not readily apparent 
that a shock forms at the cylinder surface, although this might be expected since we are 
impulsively staring the cylinder and M=0.4 is very close to the critical Mach number for 
the cylinder. From t = 2.0 to t = 10.0, in figures 3b-f, the shock formation and subsequent 
forward motion is easily seen. As the shock reaches the top of the cylinder it weakens as 
it continues to move forward until it finally disappears. The surface pressure distribution 
then tends toward the steady potential flow solution. Since the steady state solution is 
unique, it does not depend upon the time history, but the transient flow field and the work 
required to maintain the cylinder velocity during the transient phase vary greatly with the 
manner in which the steady state flow is achieved. This will be discussed in more detail 
later. 

Another interesting point to consider is that when a shock wave is present in the flow 
field, in general it will not be isentropic and hence by Crocco’s theorem (see for example 
Liepmann and Roshko [41]) vorticity may be generated. The existence of either vorticity 
or entropy in the field is of interest since the generation of these cannot be accounted 
for in an inviscid fluid without nonlinearities. Plotting the vorticity contours in figure 4 
confirms the supposition that vorticity is generated. In the figures we see that from the 
time the shock forms at t = 4.0 until the shock disappears at nearly t = 8.0, a vortex forms 
and grows in strength. After the shock disappears, the vortex follows the cylinder surface 
and subsequently convects with the fluid away from the cylinder. As the cylinder travels 
away from the vortex, the vortex diffuses as a result of its interaction with a symmetrical 
vortex of opposite sense on the other side of the symmetry plane. The vorticity magnitude 
plotted in figure 4 was calculated numerically from the dependent variables using the 
program PLOT3D [42]. 
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a) t=0.1 1 b) t=2.0 



c) t=4.0 d) t=6.0 



e) 1=8.0 f) t=1 0.0 


Figure 3. This figure is a time sequence of the pressure coefficient on the cylinder surface 
for a cylinder impulsively started to Mach 0.4. The pressure coefficient is defined as 
C p = 2(p - p 0 )/{p 0 v 2 ). 
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Figure 3. (Continued) This figure is a time sequence of the pressure coefficient on the 
cylinder surface for a cylinder impulsively started to Mach 0.4. The pressure coefficient is 
defined as C p = 2(p — p 0 )/(p 0 v 2 ). 






0 5 10 0 5 10 

e) 1-8.0 f) t-10.0 


Figure 4. This figure shows the vorticity magnitude contours, for a cylinder impulsively 
started to Mach 0.4, at various times. 
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k) t-20.0 I) t-22.0 

Figure 4. (Continued) This figure shows the vorticity magnitude contours, for a cylinder 
impulsively started to Mach 0.4, at various times. 
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Energy Partition. 

Now that we qualitatively understand the physics of the startup process, it is useful to 
discuss the problem more quantitatively in order to compare various cases. Longhorn(25/ 
has done this by examining the work required to accelerate a sphere to low Mach number 
translation using a variety of acceleration rates. If we look at the work required to acceler- 
ate the cylinder, we can develop an energy balance which describes the distribution of the 
various components of energy as a function of time. Energy flux out of some control volume 
will give us a measure of acoustic energy flux and energy flux due to convection processes. 
For the impulsively started cylinder, we will be able to distinguish between these types 
of energy by their propagation speeds, since clearly the acoustic energy will propagate at 
the fluid velocity and sound speed combined rather than the convection velocity, i.e. the 
velocity of fluid particles relative to the cylinder. 

Myers has recently developed an exact energy corollary (33/ which is well suited for 
our consideration of energy transport in flow fields which may include shock waves and 
vorticity. Myers’ result is a generalization of the concept of acoustic energy, which he 
recognized from a perturbation expansion of the general energy equation of fluid mechanics. 
Since this corollary is exact, it is uniquely appropriate for the nonlinear problems we are 
considering. For our case of perturbation about an undisturbed inviscid fluid, we can write 
Myers’ corollary as 

§ + V/ = 0. (2) 

where 

e = p(h-h 0 + — -T 0 {s-3o))-{p-j>o ) (3) 

and 

u 2 

T = p(h - h 0 + — T 0 (s — s 0 ))q. (4) 
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In these equations, h is the specific enthalpy, e + p/p, and s is the specific entropy. Recall 
that the subscript o refers to the value to of the quantity in the undisturbed medium and 
that the internal energy e and the enthalpy h are c v T and c p T, respectively. This is a 
substantially simplified version of Myers’ result since we have no mean flow, viscosity, or 
heat conduction effects. 

We can better understand this corollary if we consider each of the terms separately. 
Our starting point will be the conservative form of the total energy equation, which is 
given as 

-fif - ' + V ' ((P E +.P)l) = 0 (5) 

when pE is the total energy per unit volume, p(e + u 2 /2). This is the last equation in (1) 
if we use fixed coordinates. We need to recognize that the statement of the conservation 
of energy (5) is not unique since we can add any multiple of the continuity equation to it. 
This fact can be used to our advantage if we subtract h a times the continuity equation to 
account for energy changes in the open system which are due only to mass flux across the 
control volume boundary. Doing this, the energy equation may be written 

d u ^ 

dt + y ~ h °^ + v ‘ M* ~ h °)*) = °- ( 6 ) 

If we recall that an alternative form of the energy equation is 

Q 

Qj{ps) + V ■ {psq) = 0 ( 7 ) 

when written in terms of the specific entropy a, we can likewise account for energy changes 
due to changes in entropy by writing the energy equation as 

§i ( p I e + \ ~ h ° ~ T °( 3 ~ s °)]) + V ‘ (d* + Y~ h °~ T °( 3 ~ a °)]9) = 0- (8) 
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This is exactly Myers’ energy corollary if we add dp 0 /dt, which is identically zero. This 
energy equation will be the basis for defining the energy densities 


pu* 
tk = — 

(9) 

Ep = p(Ah — T 0 As ) — A p 

(10) 

E a — pT 0 As 

(11) 


which are the kinetic, potential, and entropy energy densities, respectively. The kinetic 
energy has its normal physical meaning and is the only component of energy in the incom- 
pressible limit. The potential energy is related to the compression of the fluid, and the 
entropy energy is energy corresponding to the increase of entropy in the fluid. We can use 
these definitions to rewrite the energy balance as 

Q 

+ £p + £ a ) + V • ((£fc + £p + e 3 )(q- v) + pq) = 0 (12) 


for a moving control volume. If we integrate over a volume V and use the divergence 
theorem, we find that 

E k +E p + E,) + F a + F a = W (13) 

when we define 



E v = 
E 9 = 


J" iv 

v 

h’ dV 

V 


(14a) 

(144) 

(14c) 
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t 

F a = J J [(ejfc + £>)(?« - Vn) + pq n \ dSdr 

t 0 g=0 


t 



t 0 g = 0 


dSdr 


w = j J 

to f = 0 


pv n dSdr. 


(1 id) 


(14e) 


(14e) 


Here E *, E p , and E 3 refer to the total amount of kinetic, potential, and entropy 
energy in the volume at a particular time. F a and F s are the total flux of acoustic energy 
and entropy energy out of the volume since time t = t 0 , and W is the total work done by 
the body on the fluid since t = t 0 . The dot over F a , F a and W in equation (13) represents 
time differentiation. The integration on / = 0 is over the body surface and on g = 0 is 
over the outer surface of the control volume V, see figure 5. Notice that we have used 
the outward normal to the body rather than to the volume in the definition of W along 
with q n = v n on / = 0. The value of F a is only truly acoustic energy flux if the control 
volume boundary is sufficiently far from the body to be relatively unaffected by the local 
aerodynamic field. We have not separated the acoustic flux into its kinetic and potential 
energy components because of the additional work term in the defining integral. This is 
normally not any problem since we are only interested in the toted acoustic energy flux. 


Time History of the Energy Transport. 

We have defined these energy quantities in order to gain an understanding of the 
transport of energy from the body surface to the near and far fields. In figures 6, 7, and 8, 
we can observe the transport of kinetic, potential, and entropy energy for the impulsively 
started Mach 0.4 translation by ex ami ning the time sequence of the appropriate energy 
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g=o 


Figure 5. The geometry of the cylinder and control volume used in the calculation of the 
toted energy in the control volume. Notice that the normal vector on / = 0 is an outward 
normal to the body and not to the volume V. 

density contours. Note that in each of these figures the energy is scaled by r to account 
for the cylindrical spreading. 

In figure 6 we clearly see that the kinetic energy is concentrated both near the cylinder 
surface and in the acoustic wavefronts. The shock is evident in the local kinetic energy 
contours as a discontinuity in the contours near the cylinder, from t = 4.0 to t = 8.0, in 
figures 6c-e. The shock becomes even more apparent when we compare figures 6c-e with 
the smooth distribution at t = 22.0 in figure 61. Notice that the local kinetic energy in 
the aerodynamic field distributes itself uniformly around the cylinder. This is not true for 
the potential energy in figure 7. Again the shock wave is evident, but now we see that the 
potential energy, which is due to compressibility, is primarily found in the wave fronts and 
near the top of the cylinder in the region of maximum fluid velocity. 

We can also see from figures 6 and 7 that more kinetic and potential energy is radiated 
in front of the cylinder than behind. This effect is at least partially numerical, due to the 
increased wave speed of the rightward moving expansion wave. We see evidence in figure 
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Figure 6. The time history of the kinetic energy field e* is shown in these figures for 
a cylinder which has been impulsively started to Mach 0.4 translation. Notice that the 
energy has been scaled here by r. 
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Figure 6. (continued) The time history of the kinetic energy field e* is shown in these figures 
for a cylinder which has been impulsively started to Mach 0.4 translation. Notice that the 
energy has been scaled here by r. 
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Figure 7. The time history of the potential energy field tp is shown in these figures for 
a cylinder which has been impulsively started to Mach 0.4 translation. Notice that the 
energy has been scaled here by r. 
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Figure 7.(continued) The time history of the potential energy field e p is shown in these 
figures for a cylinder which has been impulsively started to Mach 0.4 translation. Notice 
that the energy has been scaled here by r. 
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Figure 8. The time history of the entropy energy field e, is shown in these figures for 
a cylinder which has been impulsively started to Mach 0.4 translation. Notice that the 
energy has been scaled here by r. 
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Figure 8.(continued) The time history of the entropy energy field e s is shown in these 
figures for a cylinder which has been impulsively started to Mach 0.4 translation. Notice 
that the energy has been scaled here by r. 
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8, especially 8d, that some of the energy in the rightward moving expansion has been 
transferred to the entropy form. This should not happen for a inviscid fluid; hence it must 
in some way be related to the numerical viscosity of the discrete algorithm. Upon closer 
examination of figure 8 we find that the most intense feature roughly corresponds to the 
vortex shown in figure 3. This should have been anticipated from Crocco’s theorem. What 
is not expected in figure 8 is the elevated level of entropy energy which is clearly due to 
the passage of the acoustic waves. The entropy field has features which correspond to the 
oscillations in the density perturbation in figure 2. These figures seem to suggest that the 
numerically induced errors or damping are identified by the entropy energy term. This is 
significant since it gives a quantitative local measure of the numerical error which we can 
use to judge the acceptability of the CFD solution for acoustics. We shall consider this in 
some detail as we proceed. 

Now that we have some idea of how the components of energy are spatially and 
temporarily distributed, it is useful to consider the balance of energy for the control volume 
as a function of time. We shall consider the total kinetic, potential, and entropy energy 
in a cylindrical control volume of radius 10 around the cylinder along with the acoustic 
and entropy energy flux out of the control volume for different starting scenarios. The 
definitions for these quantities are given in equation (14), and each of the values are 
nondimensionalized by the total kinetic energy of the incompressible limit, p7ra 2 u 2 /2, which 
is the added mass times v 2 /2. In figure 9 and following, the energy components are added 
up such that the envelope of the curve represents the total work input into the fluid by the 
cylinder. Notice that the kinetic energy, potential energy and acoustic energy flux are the 
first three components. They are plotted in this order because any acoustic energy will 
be counted as kinetic and potential energy while inside the control volume and acoustic 
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energy flux as it leaves the control volume. Likewise entropy energy and entropy energy flux 
are plotted together for the same reason. Hence when the cylinder is impulsively started 
to Mach 0.4 in figure 9a for example, we see that most of the work goes initially into 
kinetic and potential energy modes and eventually a significant entropy energy component 
before t = 10.0, at which time acoustic energy begins to leave the control volume. By 
time t = 40.0 all of the acoustic energy has left and the entropy energy begins to leave. 
The entropy energy convects with the fluid velocity rather than the sound speed, thus 
accounting for the delay in leaving the control volume. We can predict the time in which 
acoustic energy and entropy energy will leave the control volume based on a local sound 
speed which is approximately unity and a convection velocity of 0.4. 

Notice that some of the energy input into the fluid by the cylinder is subsequently 
recovered as the steady state is approached. This overshoot in energy input can also be 
found in Longhorn’s calculations for the spher e[25]. We can also see that the total work 
required to reach steady state after an impulsive start, i.e. the sum of all the energy 
components at laxge t, is equally divided between the sum of kinetic and potential energy 
in the control volume and the other energy components. This result is an extension of the 
equipartition of energy found by Taylor [39], Longhom[25j, Ffowcs Williams and Lovely [38], 
and Hill/40/ since we have no low Mach number restriction and the local flow field has a 
significant potential energy component. In figure 9a, we find that the equipartition of 
energy exists even though a shock and vortices form during the transient startup process. 
These phenomenon have not been found in the previous work even though there is some 
evidence following (nonnegligible entropy energy) that a shock may be forming for Mach 
numbers as low as 0.2. 
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e) Acceleration - 2/125 


f) Acceleration = 1/125 


Figure 9. Time history of total kinetic energy E^, potential energy E p , acoustic energy 
flux Fa, entropy energy E ,, and entropy energy flux F, for a control volume with a radius 
of 10. Energy is nondimensionalized by the total incompressible kinetic energy, pira^v^ /2, 
where v is the steady state velocity Mach 0.4. 
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Figure 9 b-f shows the energy transfer for a cylinder accelerated at various constant 
acceleration rates. The fined velocity is Mach 0.4 in each case. As we compare the ac- 
celerations in figure 9, we find that the amount of work required to reach steady state 
declines significantly with reduced acceleration. As the rate of acceleration is decreased, 
the production of entropy and acoustic energies decline until in figure 9f only the work 
needed to energize the local aerodynamic field is input into the system. Thus a negligible 
amount of acoustic energy leaves the control volume. 

We find similar results for lower Mach numbers in figures 10-12. When the cylinder 
is translated at Mach 0.1 in figure 10, we see that the kinetic energy in the aerodynamic 
field is essentially that for an incompressible fluid. The compressible energy radiates as 
sound, although part of the energy we expect to radiate is contained in the entropy term. 
This is thought to be due to the numerical damping, which we will describe in more detail 
later. Also notice the oscillations in the entropy energy near time t = 10.0, in figure 
lOa-d. These oscillations are not apparent when the work is calculated directly from the 
surface pressure. The error here is attributed to the inaccuracy of the numerical method 
at this Mach number. Recall that the accuracy of the numerical scheme depends upon the 
grid resolution and time step size, neither of which vary with Mach number even though 
the low Mach number perturbations axe much smaller. Since the values in figure 10 are 
nondimensionalized by p7ro 2 u 2 /2 = .0057T, the error is more apparent. 

As we consider intermediate velocities in figures 11 and 12, we see there are no par- 
ticular surprises. The work required to bring the cylinder into steady state translation for 
each Mach n umb er and velocity is summarized in figure 13. This figure shows that more 
energy is imparted to the system for all acceleration rates when the steady state Mach 
number is raised. This additional energy goes into both the kinetic energy and the the 
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Figure 10. Time history of total kinetic energy E potential energy E p , acoustic energy 
flux F a , entropy energy E a , and entropy energy flux F t for a control volume with a radius 
of 10. Energy is nondimensionalized by the total incompressible kinetic energy, pna 2 t;^/2, 
where v is the steady state velocity Mach 0.1. 
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e) Acceleration = 2/125 


f) Acceleration = 1/125 


Figure 11. Time history of total kinetic energy E potential energy E p , acoustic energy 
flux F a , entropy energy E g , and entropy energy flux F t for a control volume with a radius 
of 10. Energy is nondimensionalized by the total incompressible kinetic energy, p-na 2 v 2 /2, 
where v is the steady state velocity Mach 0.2. 
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e) Acceleration = 2/125 


f) Acceleration - 1/125 


Figure 12. Time history of total kinetic energy E potential energy E p , acoustic energy 
flux F a , entropy energy E t , and entropy energy flux F, for a control volume with a radius 
of 10. Energy is nondimensionalized by the total incompressible kinetic energy, p7ra 2 u 2 /2, 
where v is the steady state velocity Mach 0.3. 
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Figure 13. A comparison, as a function of Mach number, of work required to reach steady 
state using various accelerations. No additional work is required to maintain the translation 
after steady state has been reached for these Mach numbers. Work is nondimensionalized 
by the total kinetic energy of the incompressible case, p7ra 2 t; 2 /2. 



potential energy components as is clearly seen in figure 14. From this figure it is evident 
that the potential energy component does not become significant until around Mach 0.2, 
hence the low Mach number solutions should be valid up to this Mach number. Again in 
figure 13, we find for all Mach numbers 0.1 to 0.4, twice as much energy is required to 
impulsively accelerate the cylinder as is required when we accelerate it very slowly, thus 
for the impulsive start an equipartition of the energy in evident. 

When the final cylinder velocity is above the critical Mach number 0.4, the character 
of the energy time history changes and no equipartition exists. Figure 15 shows the energy 
time history for Mach 0.5 translation. This case is fundamentally different from the others 
since a shock exists in the steady flow. This alters the surface pressure distribution in such 
a way that energy is continually input into the system to maintain the cylinder velocity. 
This additional energy goes into the entropy energy mode and then is convected out of the 
control volume, thus the flux of entropy continues to increase with time. The oscillations 
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Figure 14. A comparison of kinetic and potential energy in the steady state, local aero- 
dynamic field as a function of Mach number. A control volume of radius 10 is used for 
these calculations and energy is nondimensionalized by the total kinetic energy of the 
incompressible case, p7ro 2 v 2 /2. 


between the entropy energy and entropy energy flux are an indication of concentrated 
vorticity convecting out of the control volume. 

Deceleration of the Cylinder. 

Now that we have examined the energy transport from an accelerating cylinder into 
the surrounding fluid, it is informative to examine the related problem of a decelerating 
cylinder. For a cylinder impulsively brought to rest from a low Mach number translation, 
we would expect all of the energy contained in the local aerodynamic field to radiate 
as sound. For our compressible calculations, we must also entertain the possibility of 
energy going to the entropy mode. We will only consider two Mach numbers and a few 
decelerations because the results are all somewhat similar. We will plot the time histories 
of the total kinetic energy, potential energy, etc., as before, but instead of using the final 
velocity for the incompressible kinetic energy used in nondimensionalization, we shall use 
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c) Acceleration = 2/125 


Figure 15. Time history of total kinetic energy E^, potential energy E p , acoustic energy 
flux F a , entropy energy Eg, and entropy energy flux F t for a control volume with a radius 
of 10. Energy is nondimensionalized by the total incompressible kinetic energy, p7ra 2 v 2 /2, 
where v is the steady state velocity Mach 0.5. 


the initial velocity. In figures 16 and 17 we show various decelerations in which the cylinder 
is brought to rest from initial velocities of Mach 0.2 and Mach 0.4, respectively. Upon first 
glance we notice that while we have entropy energy generated in the rapid decelerations, 
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c) Deceleration = 1/25 d) Deceleration = 1/125 


Figure 16. Time history of total kinetic energy E potential energy Ep, acoustic energy 
flux F a , entropy energy E t , and entropy energy flux F t for a control volume with a radius 
of 10. Energy is nondimensionalized by the total incompressible kinetic energy, /2, 

where v = 0.2 is the velocity before the deceleration is initiated. 


there is no apparent flux of entropy energy out of the control volume. This is because the 
transport of entropy energy depends upon fluid convection, a mechanism which no longer 
exists when the cylinder is brought to rest. 

Another feature of this problem, seen for high deceleration, is the nearly linear increase 
in entropy energy before the acoustic waves leave the control volume. Since the entropy 
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Figure 17. Time history of total kinetic energy E potential energy E p , acoustic energy 
flux F a , entropy energy E t , and entropy energy flux F, for a control volume with a radius 
of 10. Energy is nondimensionalized by the total incompressible kinetic energy, pira 2 v 2 /2, 
where v = 0.4 is the velocity before the deceleration is initiated. 


energy component does not continue to grow after the acoustic waves have left the control 
volume, we may surmise that the acoustic energy is being transferred to the entropy 
term in our numerical calculations. This is true of both the kinetic and potential energy 
components and 
look back to the 


the effect is much less pronounced in the low deceleration cases. If we 
accelerating cylinder, we find evidence of the transfer of acoustic energy 
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to entropy energy in the rapid acceleration cases as well. We speculate that this is due 
in part to the inability of the discrete computational grid to resolve the high frequency 
content of the acoustic waves and therefore the numerical method damps out that part of 
the solution. Hence, the nonphysical damping arising from the discrete nature of our grid 
is seen through an increase in entropy. We will examine this idea by using different grid 
resolutions in the next section. 

A fined feature we cam understand in figures 16 and 17 is the rapid increase in potential 
energy shortly after the cylinder has been stopped. In both figure 16 a-c and 17 a-c, we can 
see that the potential energy increaise is the same as the kinetic energy decrease, aiside from 
the loss to entropy energy. We can understamd this if we consider that the fluid around 
the cylinder has inertia, and hence will interact with the cylinder as both are stopped. 
Thus, amy fluid adjacent to the cylinder will tend to stop, transferring its kinetic energy 
into potentiad energy. 

Sources of Entropy 

Before we conclude, one nagging question must still be addressed - How much of the 
entropy energy results from numerical error amd how much is real? For a cylinder which 
is rapidly stopped, it is clear that the acoustic waves transfer some of their energy into 
the entropy energy term. This is unfortunate since it implies that as an acoustic wave 
propagates, its energy will be gradually dissipated when it should continue to propagate. 
This damping of acoustic waves is evident if we examine the amount of acoustic and entropy 
energy flux out of various size control volumes. In figure 18, we find that for a control 
volume radius of 10 the acoustic flux out is greater than the entropy energy flux, while the 
opposite is true for a control volume of radius 25. This is clearly a serious problem if we 
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Figure 18. A comparison of total energy in control volumes of different radii for a cylinder 
impulsively started to Mach 0.4 translation. 


intend to use CFD methods to calculate the acoustic field directly as part of an unsteady 
calculation. 

We do have some hope that this will not be an insurmountable problem since, for lower 
acceleration rates there is very little energy in the entropy component, and hence even if 
acoustic energy i6 dissipated, it may be at negligible levels. To address the question of 
entropy energy generation, we will consider the original problem of a cylinder impulsively 
started to Mach 0.4 for the rest of our discussion. In this case some entropy energy is 
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unavoidable because of the shock wave formation, but we can also envision increases in 
entropy due to the explicit artificial dissipation added for numerical stability as well as the 
implicit damping due to the discrete nature of the numerical scheme. We shall examine 
both of these entropy generation mechanisms separately. 

Effect of Artificial Viscosity. 

The artificial viscosity we add for numerical stability, which is described in appendix 
B, seems to be an obvious possible generator of unwanted entropy energy. If we look back 
over the previous calculations, we see that the total entropy energy E a in the control volume 
is negligible when the rate of acceleration is small and it does not grow after the acoustic 
waves have left the domain. From this we should probably conclude that the second order 
dissipation is most likely to be the more important part of the artificial viscosity since 
the fourth order term does not seem to be producing any significant amount entropy after 
the gradients have left the control volume. We can test this hypothesis by using different 
values of the coefficients ^2 and £4 which control the second and fourth order dissipation 
terms, independently. In figure 19 the effect of changing the damping coefficients on the 
total potential and entropy energy is shown. We find that even if we triple k 2 or increase 
£4 by an order of magnitude, the effect on both energy components is minimal. 

Effect of Grid Resolution. 

We also expect that the grid resolution may be a strong factor since the impulsive 
accelerations and decelerations show the most significant levels of entropy energy. To study 
this effect, three additional computational grids have been generated. The grid used in the 
previous calculations and shown in figure 1 is designated here as grid 1 . Figure 20 shows 
a comparison of A r vs r for each of four grids. Grid 2 is a coarser grid in both r and 9 
directions while grids 3 and 4 both have significantly finer resolution in the r direction. 
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Figure 19. This figure shows the total potential and entropy energy in the control volume 
after a cylinder has been impulsively started into Mach 0.4 translation. The coefficients 
of the explicit artificial damping were set at: a) k 2 = .25, = .003 - the nominal values; 

b) k 2 = .50, fc 4 = .003; c) k 2 = .75, Jfe 4 = 003; d) k 2 = .25, Jb 4 = .006; e) k 2 = .25, 
jfc 4 = .020. The energy is nondimensionalized by the total steady state incompressible 
kinetic energy. 

directions while grids 3 and 4 both have significantly finer resolution in the r direction. 
Grid 3 has the finest azimuthal resolution. The constant radial spacing of the grid 4 starts 
at a much smaller r than the other grids, giving it the finest radial resolution. 

In figure 21, we compare the total potential energy and entropy energy for each of the 
four grids. Notice that there is some variation in the potential energy for both the impulsive 
case and the acceleration of 1/5 but there is a substantial variation in entropy energy. We 
can conclude that the increased entropy energy for the coarser grids corresponds directly 
to a decrease in acoustic energy flux away from the cylinder. In figures 22-24, we can 
observe the spatial distribution of kinetic, potential and entropy energies at time t = 8.0 
for each of these grids. Notice that for grid 3 the increased azimuthal resolution gives a 
smooth azimuthal variation in energy and both grids 3 and 4 have little entropy energy 
associated with the acoustic waves. These calculations support the idea that increased 
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Figure 20. This figures shows the distribution of A r for the four grids considered. The 
azimuthal resolution Ad is constant for each grid and is given as: A 6 = .0331 radians for 
grid 1; Ad = .0662 radians for grid 2; A 8 = .0165 radians for grid 3; A 8 = .0331 radians 
for grid 4. 

resolution of the computational domain enables the numerical method to better capture the 
high frequency content associated with rapid accelerations and decelerations. Hence, grid 
resolution is of primary importance in the calculation of acoustics with CFD methodology. 

Although we have refined the grid and achieved improved results, a polar grid may not 
be the best choice since the azimuthal arc length of the cells grows linearly with r. Other 
types of grids which maintain the grid cell size uniformly throughout the computational 
domain, and perhaps even unstructured grids, may be more suitable for acoustics applica- 
tions. The local aerodynamic field seems to be much less sensitive to the grid resolution 
away from the body - an observation which explains why stretched grids don’t seem to 
cause problems in aerodynamic calculations. 
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b) Acceleration 1/5 

Figure 21. A comparison of potential energy and entropy energy for various grids and 
accelerations. The final translation velocity is Mach 0.4. 

Conclusions 

The aim in part B has been twofold: 1) to understand the nonlinear, transient energy 
transfer from the surface of an accelerating cylinder to the far field; and 2) to better 
understand the capabilities and limitations of present day CFD methodology as applied to 
acoustic problems. Both of these aims have be successfully achieved through the numerical 
study of the circular cylinder. 
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Figure 22. A comparison of the kinetic energy t * distribution for various grids at t = 8.0. 
This is for the case of a cylinder impulsively started in Mach 0.4 translational motion. 


We have found that the transition from rest to a nonnegligible Mach number is some- 
what more complicated than in the low Mach number problems studied previously. When 
the cylinder is accelerated rapidly, a shock may form in the compressible fluid generating 
entropy and vortidty in the early stages of the motion. As the shock disappears, the vor- 
tidty convects away from the cylinder and the steady flow around the cylinder becomes 
essentially potential, at least for steady state velocities less than the critical Mach number 
of the cylinder. We have also shown that as the steady state Mach number increases toward 
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Figure 23. A comparison of the potential energy e p distribution for various grids at t = 8.0. 
This is for the case of a cylinder impulsively started in Mach 0.4 translational motion. 


the critical Mach number, ~ M = 0.4, the aerodynamic field experiences an increase in 
kinetic energy as well as the addition of potential (compressible) energy. These details are 
lacking in the low Mach number situation, as demonstrated by the Mach 0.1 calculations. 

It is somewhat surprising then that we are able to extend, up to the critical Mach 
number, the notion of an equipartition of energy between the local kinetic energy and the 
acoustic energy for an impulsively started body. Nontheless, our calculations show that if 
we consider the balance of energy in the local aerodynamic field, i.e. the combined kinetic 
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Figure 24. A comparison of the entropy energy t, distribution for various grids at t = 8.0. 
This is for the case of a cylinder impulsively started in Mach 0.4 translational motion. 


and potential energies, with the energy lost to the far field, i.e. the acoustic and entropy 
energy, we indeed have an equal division of the energy input into the system. Above 
the critical Mach number no such equipartition is possible because energy is continually 
transported from the body to the far field. This is accomplished through the entropy 
energy component and is due to the presence of the shock in the steady state aerodynamic 
field. 
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We have also found that separating the energy into its kinetic, potential, and entropy 
energy components is useful in understanding both the physics of the problem and the 
effect of the numerical damping. The entropy term is especially useful because it gives a 
quantitative measure that we can use to compare the effect of algorithms and grids on the 
transient solution. This use of energy is apparently new. 

In part B we have shown that for the cylinder the implicit damping in the numerical 
scheme, which is related to the resolution of the grid, is of primary importance for the 
calculation of the time dependent acoustic field. If the grid is too coarse, acoustic energy 
is transferred to entropy energy as it propagates. We have also found that the role of 
the explicitly added artificial viscosity has a small effect on the time dependent energy 
solution. 
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Summary 


The first goal of this work was to address the importance of the quadrupole source 
term in the FW-H equation. This was necessary since the acoustic analogy based theory 
has been somewhat abused in application. We concede that the quadrupole source term 
is both difficult to interpret and to evaluate. Nevertheless, we have shown in Part A that 
the quadrupole source contains fundamental components of the complete fluid mechanics 
problem, which are ignored only at the risk of error. One of the contributions of this 
research was to demonstrate that the aerodynamics problem may be solved completely only 
if the vorticity components contained in the quadrupole are retained in the formulation. 
An integral equation has been put forward which, if extended for compressible flow, could 
be useful for aeroacoustic and aeroelastic computations. With the results of Part A and 
previous investigations pointing out the importance of the quadrupole, it should be clear 
that any application of the acoustic analogy should begin with all of the source terms in 
the Ffowcs Williams and Hawkings theory. 

The second emphasis of this work was to consider the direct calculation of the acoustic 
field as part of the complete unsteady fluid mechanics problem using CFD. The recovery 
of the equipartition of energy principle for low Mach numbers gives us confidence that 
the numerical calculations are consistent with known analytical results. The finding that 
this principle can be extended to nonnegligible Mach numbers is a bonus. Now with some 
confidence, we have shown that aeroacoustic calculations can indeed be marie with CFD 
codes. Having said this, the results indicate that the acoustic field is the most susceptible 
component of the computation to numerical error. Therefore, the ability to measure the 
damping of acoustic waves is absolutely essential both to develop new algorithms and to 
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evaluate acoustic computations. Entropy energy seems to be a good quantitative measure 
for this purpose. 

This dissertation has laid essential groundwork for a new approach to the problem 
of sound generated by moving bodies. Clearly the circular cylinder solution has limited 
application in and of itself. Nonetheless, one of the main attractions of CFD technology 
is the ability to treat a wide variety of nonlinear problems with complicated geometries. 
Now with a tool to identify the numerical damping effect separately from the acoustics, it 
will be possible to develop improved numerical methods for acoustic applications. In the 
future other numerical algorithms, especially those with high order of accuracy, should be 
evaluated using model problems like the circular cylinder. New algorithms may actually 
take advantage of the entropy generation of the discrete equations by adaptively modifying 
the calculation in various spatial regions. Also the explicit damping term might be tailored 
to offset the implicit damping of the algorithm. Finally, it would be useful to separate the 
numerically generated entropy energy from that of physical origin. This might be done by 
considering the vorticity in the fluid. This new computational acoustic approach holds the 
promise of solving many problems hither to pushed aside. 
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Appendices 


Appendix A1 - Calculation of FW— H Source Terms for the Cylinder 

The incompressible form of the FW-H equation has three source terms and may be 
written 

V 2 (p'tf(/)) = -^(pv n 6{f))+-^(p'niS(f)) - ^ - (puiUjH(f)). (41.1) 

This is equation (3) of Part A. The solution may be obtained by using the two dimensional 
Green’s function for the Laplace equation in unbounded space, G(x, y) = In |x — y|/27r. In 
what follows, each source term will be evaluated analytically and by linear superposition, 
all of the contributions can be added at the end to form the final solution. The subscripts 
t, l, sind q will refer to thickness, loading, and quadrupole source terms, respectively. Note 
that since p is a constant, it will not be carried along in the intermediate calculations for 
convenience. 

Thickness Term. 

Each term of the FW-H equation may be treated independently since the wave op- 
erator is linear, therefore the thickness source will be treated first. The contribution from 
the thickness source is actually the solution to the equation 

v Vi = (.41.2) 

The solution to this may be written formally using the Green’s function as 

OO 

2*p[ = - J In I* - y \dy. (.41.3) 

— OO 
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Now since the operator V 2 commutes with d/dt it is actually beneficial to write (A1.3) as 

OO 

2 * p ' t = ~§t (^ i4 ) 

— OO 

where we define p = |x — y|. To evaluate this integral, it is appropriate for the circular 
geometry to use polar coordinates, i.e. dy = RdRdO. Now 8(f) = 6(R—a) and y = x. 0 +Rn 
if the cylinder center is labeled x 0 (t). With polar coordinates, (R, 6) are the coordinates 
of the source position y, relative to the cylinder center. Similarly, we will define the 
coordinates of the observer position x relative to the cylinder center as (r, 9). See figure 
Al-1. We can determine from the geometry that 

(? = r 2 (l — 2a cos(0 — 9) + a 2 ) (Al.5) 

where r = |x — y| and a = R/r. It is useful to recall that v n = v\ cos 9 + t »2 sin 9 and 
v = dx. 0 /dt. These definitions are used throughout the integration of each of the source 
terms. 



Figure Al-1. This figure shows the definitions of the variables used in the equations. 
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It is trivial to evaluate the R integration in equation (A1.4) due to the presence of 
the delta function. The resultant integral is 

27T 

0 f 

2tt p' t = — — / o(v i cos 9 + V 2 sin0)lnp dd (A1.6) 


which may be written as 


2ir 


, d f o 2 (v i sin 9 — V 2 cos 6) sin(0 — 9) , Q 
" — — 1 uu 


2 ' p; = I / 


r(l — 2acos(0 — 0) + a 2 ) 


(A1.7) 


after an integration by parts. Notice that now a = a/r since we are only dealing with 
the cylinder surface. Next we make a change of variables ip = 9 — 9 and recognize that 
9 < 9 < 2tt + 9 is an equally valid range of integration in (Al.7). After some algebra we 
can write 


2 tt 


2-npt 


f d [a 2 (v„ sin ip — cos ip) sin ip 


.*[ 

dt J 


r(l — 2a cos ip + a 2 ) 


dip 


(A1.8) 


where v ^ = v\ cos 9 + v 2 sin 9 and vj = —v\ sin 9 + v 2 cos 9. In Gradshteyn and Ryzhik [43], 
we find 

2ir 


J~- 


sin(nip) 


2a cos ip + a 2 

cos(nip) 

2a cos ip + a 


dip = 0 


2ir 

cos^ny/; J, _ *71 « 2 

1- 


27ra n z 

2 dip = - — -5- ; at < 1. 


1 — a 2 


(A1.9) 


(A1.10) 


Clearly a < 1 will always be true for this case since a = a/r and a is the minimum radius 
of interest. Substituting these results into (A1.8) gives 


2 d ,v • n, 


27r P ; = a 2 7r^(— 


(Al.ll) 
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with v • n = Vfi and ii = (x — x 0 )/|x — x G |. Now using the relations 


d 1 v • n 
dt r r 2 


and 


d_ 

dt 


1 (v • n) . 
r 


v . _ \ X 

_( n ) = --v+ - " n 


(A1.12) 


(A1.13) 


the solution for the thickness contribution to the pressure perturbation is found to be 

p ' (x ' t) = 2 ~ v '‘ ) *TTt ' f } < A114 ) 

when written in terms of the usual variables and including the fluid density p. 

Loading Term. 

The loading term may be determined in a manner similar to the thickness term, 
using the above definitions. After the initial manipulations and completion of the radial 
integration, the loading term is found to be 

2tt 

d f 

2irp'i = — — / ap| njlnp dd (Al.15) 

OXi J r—a 


where 


, „ 2 3v 2 dv „ K K 2 

p| = 2v n - IT + °^7 • n - — y t - 


r—a 


(A1.16) 


2 ’ ~dt ~ 7T a 1 8tt 2 o 2 ' 

Now we can take the spatial derivatives inside the integrand where they act only on p. 
Noting that 

dp ( xi - x 0i - an { ) 


dx{ 


and 


dp r(cos(0 — 9) — a) 

P 


ni dxi 


(A1.17) 


(A1.18) 
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we can write equation (.41.16) as 


2 ir . 


2np'i = a J 


( 2 v j + ~ I * 2 ~ lk?t ~ l£?~)( cos ( g ~ j) ~ a ) 

1 — 2acos(0 — 0) + a 2 




d0. (41.19) 


Using some trigonometry and making a change of variables puts this equation in the form 

2ir 


2n p'i = a J 


o 

2ir 


(fol ~ v i 2 ) cog ( 2 ^) + 2 ( a ir • ” - £**) cos V 1 ) ( cos - a ) 

1 — 2a cos ip + a 2 


dif) 


(v h v- t sin(2V>) + (»!*+ s ‘ n V’-” 2 - w) (cos it - «) , „ A> 

+ ° ' l-2acos,t+a* ( ’ 


°/ 


This is now in a form in which we can apply (.41.9) and (4.1.10) quite easily. We find that 
the second integral is zero, therefore the solution to loading source term contribution to 
the pressure perturbation is 


.. . p f a 2 . 2 2x a 2 dv „ Kvt 

rf(=c,0 = |{^(»S-«?) + T ^ r - — 


(41.21) 


when written in terms of the usual variables and including the fluid density p. 


Quadrupole Term. 

The quadrupole source term is slighty more complex since the radial integration is no 
longer trivial. The formed solution for the pressure due to the quadrupole is written as 
2tt oo 

2np'q = — f j — — (u{UjH(R — a)) In ^r 2 — 2rR cos(0 — 0) + ~R? RdRdd. (41.22) 
J J oxiOxj 
0 0 

One way to evaluate this integral is to recognize that 


+ + Ml-**) 
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and split the problem into three integrals 

2tt oo 


h = ~ J J ln \/ r2 -2ri2cos(0-0) + .R2 RdRdd (A1.24) 

0 a J 


27T oo 


I 2 = -^~ j J uiu n 6(R - a) ln yjr 2 - 2rflcos(0 - J) + ~R 2 RdRdd (A1 .25) 


0 0 
27T oo 


J 3 = - J J -^-ni6{R-a)ln y/7* - 2r.Rcos(0 - 0) + i? 2 RdRdO. (A1.26) 
0 0 J 


The first integral, Jj, requires more attention than the other two since we must be 
careful about the radial integration, so we will leave it until last. The second integral, 1 2 , 
is in a form which is very similar to the loading term. Here we need U{U n evaluated on the 
surface t = a. From the velocity potential we know that 

a 2 

Ui = - 5 -(v„n t - - v t ti) - - — ti (41.27) 

r* Zirr 

and hence 

K 

Ui tin = ((vnTlj - v t ti) - — (41-28) 

on the surface. Here n t - and <,• are the components of the unit normal an tangent vectors 
on the body surface. Now if we integrate radially and bring the spatial derivative in the 
integral, I 2 becomes 


h = 


f (cos(fl - 0) - a) + ( VnVt - Kvn) sin(0 - 0) 

J 1 — 2 a cos (9 — 6) + a 2 


{A 1.29) 


As before, we need to change variables, ip = 8 — 6, and use the identities v n — v „ cos ip + 
v- t simp and vt — —v^ sin ip + Vj cos ip, along with some trigonometry, to get an integral of 
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the form of (41.9) and (41.10). Skipping the details since we have already done most of 
them for the loading term, we find that 


r K 

h = - T n 


The integral J 3 is slightly different in that we will need to know that 

duiUj du { 

— ni = —u jni 

for an incompressible fluid. Using the velocity potential, we can write 1 3 as 

K 2 3 K 


(41.30) 


(41.31) 


2 tt 


r 2 K* 3 A - N , 

I * = J( 2v + i® + s”‘) lD '* 


(41.32) 


The third term in this integral is of the same form as the integral in (41.7) and may be 
evaluated immediately. The integration of the first and second terms is divergent, but it 
will be seen that these cancel exactly with part of the integral I\. For this reason we will 
write the solution to J 3 as 

2ir 

K 2 

(41.33) 


r f, n 2 Xl 3 K 

,3 =J ^ + 4^ )lnpde ~2T Vt - 


Now to solve the final integral, I \ , we first make use of the fact that for incompressible 
flow the relationship 


d 2 UiUj du{ duj 8 a 4 2 4A K 2 

~ = — 5 -v H r v t + 


dxidxj du j dx{ r 6 7rr 5 27r 2 r 4 

is true for the circular cylinder. Thus the integral to be evaluated is 


( 41 . 34 ) 


27T 00 


-//(£• 

0 a 


8a 4 •> 4 K 


+ — rv t + 


K 2 


ttA 4 ‘ ' 2tt2A3 


^ In y/7 2 -2rflcos(0 -0) + A 2 dRd9. (41.35) 
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If we integrate by parts this integral becomes 

= + + i npd e 


0 

2n oo 


+ 


f f( 2a \l I iKa \ _ I (cos(» -»)-«) 

JJkR* + 3n& ‘ + 4ir2«Vr(l-2acos(«— #) + a 2) rffi ( 1 

0 a v 

with a again defined as a = i?/r. If we call the first integral in (.41.36) I\ a , we can 
immediately write down the solution 

2ir 


=-/( 


2 X, ^ 4tf 

2v + toV )ll, ' ,<W + 37'" 


(41.37) 


by looking at 1 3. The second integral in (41.36), In,, can be solved by reversing the order 
of integration and breaking the integration in two parts to make use of the formulas from 
Gradshteyn and Ryzhik, (41.9), (41. 10), and 

2ir 


/ 


cos (nip) 


1 — 2a cos ip + a 2 


dip = 


27T 


(1 — a 2 )a n 


; a 1 > 1. 


(41.38) 


When this written out we find that 
r 00 2tt 


T / f /Y2a 4 o 4Ka 2 if 2 \ (cos(0 - 0) - a) 


cos (0 — 0) + a 2 ) 


dd. (41.39) 


Notice that for the first integral, f* dR, the range of integration restricts a to a < 1 while 
in the second integral, f£° dR, a > 1 is true. For a < 1 the 0 integration has been done 
previously in (41.19)-(41.21), with different coefficients. For a > 1 we can write the 0 
integration as 

r /Y 2 ° 4 2 K 2 4tfa 2 • , # x\ (< 

I= J\-w v + -^ + ^- VASm * +viaxi,) )^ 


COS 


Ip - a) 


2a cos ip + a 2 ) 


dip. (41.40) 
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Upon application of (A1.9), (41.10), and (41.38) to this we can write In, as 

r } AKa 2 Jn 7 n /2a 4 , K 2 . 4Ka 2 r , 

hh = J MP v -‘ dR ~ J + + -Wr v i dR - ( ' 4L41) 


It is straightforward to complete the integration in R and we find that 
T 2 f/ n 2 K 2 v, „ 7 TO 4 2 K 2 ,Ka 2 2K. 

I ' = -J( 2v + W )ln '’ d> --^ V 


(41.42) 


Combining I\, I 2 , and I 3 gives the full quadrupole contribution of the pressure perturbation 


as 


Ka 2 


K 2 


Pff(*.0 = - 3 ( 74 * + ~^T Vt + 4 ^ 2 ) 

when written in terms of the usual variables and including the fluid density p. 


(41.43) 


Appendix A 2 - Derivation of New Quadrupole Formulation 

The quadrupole source in the FW-H equation can be arranged into various forms 
which can be useful in the interpretation of the role of this source term. We can derive 
equation (17) by first expanding the quadrupole into the form 


= V ' {(p«V-u + u- V(pu))H(f) +put>„5(/)j (42.1) 


d 2 (pujU jH(f)) 

dx 


where we have used Vff (/) = n£(/) and u n = v n on / = 0. Now if we use the relation 

,1 


u • Vu = V i-u 2 4- £ x u, 


(42.2) 


where £ = V x u, and the fact that for an incompressible flow, V • u = 0 and p constant, 
we can write the quadrupole 


;"^ (/)) = V ■ j(vi*m 2 + p(xu)H(f)+ pv„uii(/) J . (A2.3) 


dx 
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Now if we take the Heaviside function inside the V^pu 2 term we get 
82( a’x’af (/>> = V 2 (^ 2 Jf(/)) + pV ■ (C x ug(/)) + V ■ ((pv n u - \pu 2 n)S(f)) (A2A) 

which is the desired result. Notice that the vector quantities in the surface term may also 
be written as 

u 2 ^ v 2 — u? « u 2 „ 

v n u-yn= - — L n + v n u t t = — m (-42.5) 

where m is a unit vector. 

Now we can also rearrange the thickness source term by some simple manipulations. 
The thickness source term can be expanded as 

I (p°«i(f)) = ■ ft. *(/) + pv • a(fi ^ (/)) . (>12.6) 

If we recognize that 

^ = -V/v (-42.7) 

and V/ = n, then with this and some algebra we can also show that 

^(ft «(/)) = -V(«n«(/)) (-42.8) 

which allows us to write equation (.A2.6) as 

a(yv " <(/)) = a - py • V(»ni(/)). (>12.9) 

Cft Cft 

Now since v = v„ + w x rj we can show that V • v = 0. Here v 0 is the velocity of the center 
of rotation, rj is a position vector from the center of rotation, and u> is the rotation vector. 
This result is clear when we write V -v = T] - V x u> — u; • V xrj because Vxw and V x rj 
axe both zero. Finally we may write the incompressible version of the thickness source as 

4 (pv n eu)) = P% ■ ft - v • (p»„v<(/)) . (-42.10) 
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This form is useful because it can be combined with the new quadrupole description to 
write the incompressible FW-H equation as 

V 2 ((p' + (/)) =V ■ {((p' + jpu 2 )A - pv „( ut - «()£)*(/)} - p-^ ■ nS(f ) 

- f V-(((xuW)). U2.ll) 

Dowling [44] has suggested an alternate derivation of equation (A2.ll) once it is 
recognized that the quatitity p' + jpu 2 is the quantity of interest. It starts with the 
conservation of momentum 

p(^ + u-Vu) = -Vp = -Vp' (.42.12) 

and the relation (A2.2). Now if B = p' + pu 2 / 2 we can write 

H{f)VB = -p(^ + C x u)H(f) (A2.13) 

which is valid in the entire space. Next if we replace H(f)VB with V{BH(f)} — BiiS(f) 
and take the divergence of equation (A2.12) we have 

V 2 (£ff(/)) = V • (BnS(f)) - p ^ ■ n 1(f) - pV ■ ((C x u )H(f)) (A2.14) 

since V • u = 0. Comping this equation with (A2.ll) shows that 

|£ • n 6(f) = ^ • nS(f) + V ■ (v n (u t - v,)U(f)). (42.15) 
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Appendix A3 - Comparison of Acoustic Theory to Classical Aerodynamics 

In this appendix it is our goal to show that equation (20) can be reduced to classical 
thin airfoil theory for both the thickness and lifting problems. To do this we need to first 
start with the differential equation (20) 

V 2 (BB(/)) = V ■ {(Bn - pv n [u, - v, ){)«(/)} (43.1) 

which has been simplified by the additional assumptions of irrotationality and steady body 
motion. If we recognize that 

v n (ut - vt) = v n ut - vtu n = v • (u x k) (A3. 2) 

we can write the integral representation of the solution as 

1 / Bcosd — pv • (u x 

B= 2 7 / 7 

/= 0 

by using the Green’s function technique. In this expression we have used cos 6 = r • n, 
sin# = r • t, and r = |x — y|. This equation is valid for a reference frame in which the 
body is moving whereas in aerodynamics the body is stationary and the fluid is moving. 
This difference in point of view is the motivation for comparing the Bernoulli equation for 
the two reference frames, where we find 

r' = + = -IpJ-fPA. (43.4) 

Here the first right hand side is the normal, body fixed result from the aerodynamics 
and the second right hand side is the fluid fixed, acoustics result. In both cases u is the 
magnitude of the velocity perturbation of the fluid. From this we can infer that 

B = p' + = p\ ■ u (A3. 5) 


k) sin 6 


dS 


(A 3.3) 
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since the constant body velocity v = — U<x,. Combining this relation with equation (A3. 3) 
allows us to write the vector equation 


1 f u cos 6 — (u x k) sin 9 

u = — / i dS + w. 

27 r y r 

/= 0 


(A3. 6) 


Here a; is a vector field everywhere orthogonal to the body velocity v. The vector field 
u> must be irrotational, Vxw = 0, and a solution to the Laplace equation, V 2 w = 0, for 
u to be a solution. This combination of conditions together with the orthogonality to v 
requires that u> is of the form 

u ) = a + bq (A3. 7) 


where rj is a distance in the coordinate direction orthogonal to v and a and 6 are constants. 
Further, a and b must both be zero since the velocity perturbation is required to approach 
zero at infinity. Therefore we can write equation (A3.6) as 


u 


1 t u cos 6 — (u x k) sin Q 

2i x J T 

f = 0 


(A3. 8) 


without ambiguity. This equation will be the starting point of our comparison with thin 
airfoil theory for the separate cases of a symmetric nonlifting airfoil and a thin, cambered 
airfoil at angle of attack. 
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Slender Body Theory. 

In the first case consider a symmetric slender two dimensional airfoil which is described 
by the profile 

y = ±F(x), 0 < x < 1 ( A2.9 ) 

and is moving along the x axis with velocity U 0 0 . If a continuous distribution of sources 
and sinks is distributed along the axis, it is well known (see Yih [32] for example), that the 
velocity potential may be written 

0 s (43.10) 

where /(£) is the source strength at x = £. The linearized boundary condition on the 
surface of the body gives 

v = ±Uoo^ (A3.11) 

ax 

where the positive and negative signs are for the upper and lower surfaces respectively. 
Upon differentiation of (A3. 10) with respect to x and y, we find that the fluid velocities 
are given as 


u= _l r m—o M 

Wo (x-0 2 + y 2 ^ 

{A 3.12) 

v= jl [' m * 

2?r Jq (x-tf+y* 

(A3.13) 


For small positive y we can evaluate an approximation to the second integral and after 
using (A3. 11) it follows that 

/(*)=2Uoo^H M3.14) 
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While this is all we need for our comparison, usually the first velocity component is given 
as the finite part of 

i r 1 m 


u = ~ f 

2tt Jq x - ( 


(A 3.15) 

jo * — ^ 

for small y. 

For the acoustic equation (A3. 8), we need the surface normal n a (—dF/dx, ±1) for 
the upper and lower surface. We also need 

(x - £)ni + yn 2 


r • n = 


f • t = 


~(x - 0n 2 + yni 


(A3. 16) 
(A3. 17) 


to be able to rewrite equation (A3. 8) as 

1 f 1 u( ~(s ~ + y) ~ (u x &)(-(» - 0 ~ 

(x - t) 2 + y 2 

1 f 1 u(-(x - Qf - y) - (u x k)((s - Q - yg) 


u= a 


dt 


In Jo 


+ -' ' ' "" (s-0 2 + /' - 1 (A318) 

Here the first integral refers to the upper surface and the second integral to the lower 
surface. Since the airfoil is symmetric, we know that u\ upper = u \lower an d v\ U p per = 
~ v \lower- Also we can write u x k = (v, — u). Combining this information and dropping 
second order terms it follows that 


U= 1 f m—q M 

2n Jo (x - £) 2 + y 2 

(A 3,19) 

V f l 2v 

2tt Jq (x - £) 2 + y 2 

(A3. 20) 


which is identical to the slender body theory result (A3. 12) and ( A3.13) when the linearized 
boundary condition (A3. 11) is applied in terms of the upper surface. This result shows 
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the equivalence of the acoustic equation to the two dimensioned slender body theory after 
these approximations are made. 


Munk’s Vortex Sheet Theory. 

If we have a thin airfoil which is producing lift, we can model the airfoil by its mean 
line and consider it to be a vortex sheet on the x axis. If we do this, we can write the 
velocity potential 

(43.21) 




where f{£) is the vortex sheet strength at x = £. The tangential velocity on the airfoil is 
approximately 

d<f> y f 1 m 


u _ cty _ y_ r 

U dx 27r J 0 


(* - 0 2 + y 2 

which is similar to (A3. 13) and from this we can deduce 

/(*) 


dt 


(4 3.22) 


u\ = = -u\ 

upper 2 lower 


(A 3.23) 


and hence the vortex sheet strength is f(x) = [tt], the discontinuity in u. To find the value 
of /(£), we note that the hneax boundary condition gives 


TT ( dym N 

v = Uc ° ( ~dT ~ a) 


(A 3.24) 


with y m defined as the distance to the mean line from the chord fine and a is the angle of 
attack. We cam also find that 


„ = d ± = L f l JMzzJLjt 

dy 2jt 7 0 (i - 0 2 + y 2 

Therefore if we consider small y, we can write (A3.25) ats 

v = r t 

27 t Jo x-i 


(43.25) 


(43.26) 
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which is an integral equation in terms of /(£) when we consider it together with (A3. 24). 
We must solve this integral equation with the Kutta condition /(l) = 0 to ensure there is 
no infinite velocity at the trailing edge. The solution of (A3. 26) is then known to be 




{ v (0 
1 - t ~ *) 


di. 


(A3. 27) 


This brief outline follows that of Yih [32] closely. 

Now for the acoustic equation we first write the normal vector as n « (— dy m /dx+a , 1) 
along with 


r . (*-0(-% + a) + » 

r 

. 2 ~( x ~ 0 + + <*) 

r . t = 


(A3. 28) 
(A3. 29) 


for the upper surface. Also for this thin airfoil approximation, n| upper = — Slower i which 
enables us to write (A3. 8) as 

u = 2. r 1 M((» - *)(-%■ + «)+») + M x £(-(* - j) + »(- % + °)) K (A3 30) 

27 r Jq r 2 


if the discontinuity in u is written [u]. Keeping only the first order terms, it follows that 


. y r M r 

2ir Jo (* - 0 2 + y 2 

(A3. 31) 

V = 1 /' [“](*-{) ^ 
2 tt Jo ( X - O 2 + y 2 

(A3. 32) 


where we have used the fact that [v] « 0. Notice (A3.31) corresponds to (A3. 22) and 
(A3. 32) corresponds to (A3. 25). Thus we have shown that we can reduce the acoustic 
equation to Munk’s vortex sheet theory as well as slender body theory. 
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Appendix B - Description of the Euler Solver 

The purpose of this appendix is to give a detailed description of the numerical proce- 
dure used in the computations included in Part B. Some description of the methodology is 
given, but most of the concepts follow Jameson’s work. References [36-39] give a somewhat 
broader description of the ideas used in the code developed for this work, while variations 
from normal practice are pointed out here. 

Euler Equations for a Moving Grid. 

The Euler equations may be written in vectorial form as 

^ + VF = 0 (Bl.l) 


where the dependent variable vector is 


U = 


r P ' 

pu 

< > 

pv 

. pE > 


(£ 1 . 2 ) 


and the flux matrix is 


F = Vq + p < ) >• 

l q , 


(B 1-3) 


In these equations, p is the fluid density, pu,pv are the components of the momentum in 
the xi and directions, pE is the total energy ( E = c v T + u 2 /2), while the fluid velocity 
q = (u,u) and i, j are the unit vectors in the x\ and *2 coordinate directions. This form 
of the equations is known as the conservation form and is the correct form to accurately 
calculate the location and strength of shock waves numerically, if they exist in the flow 
field. 
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If we wish to consider a moving control volume, it is useful to make a change of 


coordinates 


X(<) = x + J v{r)di 


(51.4) 


where v(t) is the velocity the moving coordinate system. When we substitute the moving 
coordinates X in the Euler equations (Bl.l) we find that 


au dU „ TT 

~w ~ nr + v (*) ■ 

<n dt x(t) 


(51.5) 


Vx • F = V x • F. 


(51.6) 


Thus writing the Euler equations in terms of the moving coordinates gives us 


+ v x (F-Utr)=o. 


(51.7) 


We can arrive at the integral form of this result if we use generalized functions and extend 
the Euler equations Bl.l outside of the control volume to the entire 2-D space. If we 
consider a control volume V where g > 0, bounded by g = 0, the extended equations 


would be 


£(Uff( S )) + V - (F H(g)) - (u| + F ■ Vg)S(g) = 0. 


(« 1 - 8 ) 


We will define g such that V <7 is a unit normal vector to the surface g = 0, pointing into 


the control volume. Recall that 


^ + v-Vg = 0 


(B 19) 


for a surface with velocity v. Then if we integrate over the entire space, we find 


oo 

lLJudy+ J W (FH{g))dy+ J (F - Uu) • hdS = 0 (51.10)- 


98 



where n = -Vg, the outward unit normal and v is the velocity of the control volume 
surface. The second integral in B1.10 is identically zero, which we can see if we use the 
divergence thereom and realize that the Heaviside function is zero on the integration surface 
since g < 0 outside the control volume. Hence we may write Euler equations for a moving 
control volume in integral form as 

^ + /(**- Uv)-ndS = 0. (B 1.11) 

V g=0 

If we change the definition of F to include the Uu term in either Bl.7orBl.ll, the form 
of the equations is identical to the one we started with. Notice that in any case the fluid 
velocities are defined in the reference frame fixed to the undisturbed medium. 

The Finite Volume Approach. 

Starting with the last equation and changing F as indicated, we have 

^ J UdV + J F • dS = 0 (£1.12) 

V s 

where S is the boundary of the control volume. For a finite volume method we assume that 
we can break up the domain of interest into small volumes AV and then if we interpret U 
to be the value of the dependent variables for the volume, in the sense of the mean value 
theorem, we can write a system of ordinary differential equations for each AV. This idea 
leads to the semi-discrete version of the conservation law 

^(UAF) + ^ F • AS = 0. (£1.13) 

faces 

The surface integral has been written here as the sum of the flux through each face of AV , 
with AS an outward normal vector to the face with magnitude AS. Here AS is the surface 
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area of the particular face rather than of the complete volume AV. For two dimensions 
the volume AV should be interpreted as the area of the cell and the surface area A 5 as 
the length of an edge of the cell. The value of the dependent variables [p, pu, pv, pE]^ at 
the cell face is approximated as the average value of the two adjacent cells. The pressure 
at the cell face is directly calculated from the average value of the dependent variables. 
This scheme reduces to a 0(Ax 2 ) accurate central difference scheme on a Cartesian grid 
and has the property that uniform flow is an exact solution to the difference equations. 

One advantage of the finite volume method is that curvilinear coordinates can be used 
to develop the computational mesh, and therefore AV and AS, while regular Cartesian 
coordinates of the fluid velocity, (u,v), are used in the specification of the dependent 
variables. Hence no transformation is needed and the form remains the same for various 
types of grids. The finite volume method is also readily applicable to triangular grids and 
unstructured grids. In this work only polar grids with constant azimuthal spacing are used. 

Dissipative Terms. 

The scheme presented thus far permits high frequency oscillations between odd and 
even mesh points which may cause the solution to become unstable. Dissipative terms 
must be added to suppress these spurious oscillations and prevent nonphysical wiggles 
near shock waves. The semi-discrete form of the equations can be written 

^(Ui,AVy) + Qij - Dy = 0. (£1.14) 

where Qy is the convective flux balance vector for the i/th grid cell and Dy is the added 
dissipative term for the same cell. No summation of indices is implied in this equation. 
Various forms of D jy axe possible, each changing the character of the difference equations 
as Jameson and Lax/38/ have shown. 
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For the present work, Jameson’s combination of second and fourth order, adaptive 
dissipation has been used [37]. This form of dissipative terms is attractive because it is 
computationally efficient, robust, and in wide use. The dissipation Djy is written in this 
form as 

Dij = Afdi + Afdj (£1.15) 

where A + is the forward difference operator, A f di = \ — d{, and d t - is 

d, = ( AVi ^ Vi ~ -)(e (2) Ar u <j - £ (4 >(Ar) 3 U ii ). (31.16) 

Here A - is the backward difference operator, A“U {j = U ij — Uj_i j, and (A“) 3 means 
three applications of A ~ . The coefficients and are defined as 

= k 2 max(i/j_ 2 , 1'i-i) v i+\) (£1.17) 

and 

g( 4 ) = max(0, k\ — e^ 2 )). (£1.18) 

The pressure switch V{, which is designed to detect the presence of shocks, is defined as 

_ Pi+lj ~ ^Pij Pi— lj 

Pi+ij + 2pij + Pi-lj 

and the constants Jb 2 and k^ have been used with nominal values of k 2 = -25 and = .003. 
The value of d j is found in an analogous manner. From these definitions we see that the 
constant controls the fourth order background dissipation and k 2 controls the second 
order dissipation around pressure gradients. The effect of the second order term is to make 
the numerical scheme locally 0( Az) accurate and it is clear that the pressure sensor v , 
will interpret any high frequency pressure variation as a shock and try to dampen it. The 
effect of the coefficients of k 2 and on the solution is shown in Part B. 


(£1.19) 
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Multistage Time Stepping. 

If we keep the cell volume A V{j time independent, we can write the semi-discrete 
system of equations for each cell as 

3 (Ui ' ) + A^J (Q «~ Di ' ) = 0 - (B1 ' 20) 

Recall that Q is the convective operator and D is the dissipative operator. Multistage 
time stepping schemes are modeled after the Runge-Kutta procedures used for ordinary 
differential equations. By adding more stages we can either improve the time accuracy or 
increase the allowable time step size while maintaining stability. A five stage scheme of 
the form 


U<!> 

u< 3 > 

u< 4 > 

*J 


At 


= u 5- a iAi7r(Q5- D S) 


A Vij 

- U” - ao-— — - D^) 

XJ 2 A X} 

= m - a 3 -^-(Q( 2 > - D^) 

3 3 A Vij K ^ tJ XJ ' 

-U?_a 4 — 

- u v 4 A Vij x3 xj } 


At 


U ?. +1 = U£ - -^-(q( 4) - D^) 

V7 AVii ^ 13 13 


(51. 21a) 
(51.216) 
(51.21c) 
(51.21d) 
(51.21e) 


has been used for all the calculations in this work. Notice that only two evaluations of 
the dissipation are made in order to reduce the number of calculations. The coefficients 
ai = 1/4, at 2 = 1/6, a 3 = 3/8, and <*4 = 1/2 have been used. This time discretization 
is 0(A< 2 ) accurate. Jameson’s papers [36,37,39] are good references for more detail on the 
multistage time stepping approach. 
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Boundary Conditions. 

Numerical boundary conditions are needed at all of the computational boundaries to 
complete the system of equations. For exterior flow problems artificial boundaries are nec- 
essary to produce a bounded domain. Improper treatment of these boundary conditions 
can lead to serious errors and perhaps instabilities, thus for this work the boundary condi- 
tions recommended by Schmidt and Jameson have been used [37]. Many other formulations 
of the numerical boundary conditions exist, but it was not within the scope of the present 
work to consider a wide variety of boundary condition formulations. This section follows 
Schmidt and Jameson and is given for completeness. 

For the outer boundaries in subsonic flow, there is one outgoing characteristic at the 
inflow boundaries and three outgoing characteristics at the outflow boundaries. According 
to the theory of Kreiss [45], three conditions need to be specified at the inflow boundary 
and one at the outflow boundary. The remaining conditions required by the numerical 
method are determined by extrapolating the outgoing characteristic variables from the 
interior solution. Thus at the inflow boundaries we have used the linear characteristic 
relations 


p - clp = p 0 - c 2 0 po 

(B 1.22a) 

s 

II 

s 

0 

(£1.226) 

P — Po^oOn — Pc ~ Po c o<ln e 

(£1.22c) 

P + poCoQn = Pe + PoCoQnt 

(Bl.22d) 


where the subscripts o and e refer to the value in the undisturbed medium and the value 
extrapolated from the interior of the computational domain, respectively. Here q n and 
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qt are the normal and tangential components of velocity. When combined these give the 
result 


P=\(Pe+Po + PoC 0 (qn e ~ 9n 0 )) 

n -n 4 _P-Po 

9n ~ 9n 0 H 

Po c o 


(B1.23a) 

(B1.236) 


and the density can be determined from (B1.22a). 

. For the outflow boundary, a nonreflecting condition which should eliminate incoming 
waves is 

"q^JP ~ PoCo9n) — 0. (SI. 24) 


This condition has been used together with the extrapolation of the density and momentum 
from the interior to minimize non-physical reflections from the outflow boundary. This 
boundary condition does in fact allow small reflections in actual application, which can 
lead to instabilities, however for the calculations considered here the time required for 
acoustic waves to reach the boundary is easily determined. Thus by moving the outflow 
boundary sufficiently far away from the body and limiting the computational time, any 
problems with the outflow boundary conditions can be minimized. 

For the surface boundary conditions, no flux through the surface is specified, i.e. 


FA S = p 


0 

ni 

V n 


AS 


{B 1.25) 


since q n = v n on the surface of the body. The pressure at the surface is taken to be the 
linear extrapolation of the pressure from the interior. Linear extrapolation is used because 
of its simplicity and because little difference was seen between the results obtained with 
more sophisticated surface boundary conditions. The success of linear extrapolation is 
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probably related to the fact that the computational grid used in the present calculations 
is very fine near the body. 

The solution is assumed symmetric about the x axis in the computations presented. 
Calculations were also made with a full computational grid with no symmetry assumption 
for the circular cylinder, but the results were identical. The boundary conditions used in 
this work should be considered representative of those in common use even though more 
advanced boundary condition formulations are under development. 
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